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Preface 


The  ALR-69  radar  warning  receiver  on  the  F-16  has 
been  hampered  by  the  presence  of  secondary  lobes  vdiich 
cause  errors  in  the  direction  finding  routine.  This  study 
examined  this  problem  in  terms  of  cylindrical  wave  scatter¬ 
ing  from  a  large  dielectric  shell.  The  shell  was  modeled 
as  an  ellipse  since  this  was  a  conic  shape  that  was  close 
to  the  actual  shape  of  the  radome.  The  intent  of  this 
study  was  to  examine  the  problem  and  see  if  the  side  lobes 
are  caused  by  scattering. 

The  solution  method  was  to  use  a  Galerkin  applica¬ 
tion  of  the  method  of  moments.  The  electric  field  expan¬ 
sion  fvinction  was  the  piecewise  sinusoidal  basis  function. 
The  integral  equation  developed  was  Richmond's  Integral 
Equation  which  is  valid  for  any  dielectric  cylindrical 
shell.  The  shape  of  the  object  resulted  in  the  use  of 
elliptic  coordinates  vdiich  is  not  one  of  the  standard 
orthogonal  coordinate  systems.  The  combination  of  ellip¬ 
tic  coordinates  and  Hankel  functions  made  the  integration 
nonexistent  in  closed  form.  The  resultant  numerical  inte¬ 
gration  took  a  great  deal  of  computational  time.  This  time 
problem  was  further  compounded  by  the  presence  of  a  line 
singularity  involving  the  Hankel  function.  Considerable 


discussion  is  given  the  subject  from  both  the  mathematical 
aspect  and  the  programming  aspect. 

This  thesis  did  not  generate  any  far  field  plots 
of  the  electric  field  since  the  program  used  to  calculate 
the  reaction  matrix  elements  produced  erroneous  data.  The 
exact  reason  is  unknown.  The  theory  employed  by  this  work 
is  not  new,  it  is  only  a  different  application.  The  sinqple 
fact  that  it  took  over  1300  CPU  seconds  to  fill  a  33  x  33 
matrix  and  a  33  element  vector  on  a  machine  as  fast  as  the 
CDC  Cyber  175  indicates  an  impractability  of  the  method 
used. 

A  special  note  of  thanks  to  my  sponsors,  Mr.  William 
Kent  and  First  Lieutenant  Robert  Schneider,  ASD/ENAMA. 

The  amount  of  help  given  by  providing  me  with  an  ASD/EN 
problem  number,  and  account  for  the  cyber  is  immeasurable.. 
Thanks  are  due  to  my  advisor.  Captain  Thomas  W.  Johnson, 
who  was  personally  excited  and  motivated  by  the  research. 
Finally,  a  special  note  of  gratitude  to  Mary  Browning, 

Linda  Stoddart,  and  Veleta  Kendall,  AFIT/LDE.  These  people 
found  information  from  the  most  unusual  sources  possible, 
euid  were  a  real  help  in  getting  this  project  anyvdiere. 
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Abstract 


This  thesis  examines  the  scattering  of  cylindrical 
waves  by  large  dielectric  scatterers  of  elliptic  cross 
section.  The  solution  method  was  the  method  of  moments 


using  a  Galerkin  approach.  Sinusoidal  basis  and  testing 
functions  were  used  resulting  in  a  higher  convergence  rate 
The  higher  rate  of  convergence  made  it  possible  for  the 
program  to  run  on  the  Aeronautical  Systems  Division's 


CYBER  computers  without  any  special  storage  methods. 

The  program  thus  developed  required  very^,J^ge 


run  times.  This  makes  the  program  impractical  for  scatter 
ers  of  size  greater  than  one  wavelength.  -^This  report 


includes  discussion  on  moment  methods,  solution  of  inte¬ 


gral  equations,  and  the  relationship  between  the  electric 
field  and  the  source  region  or  self  cell  singularity. 
Since  the  program  produced  unacceptable  run  times,  no 


results  are  contained  herein.  The  importance  of  this  work 
is  the  evaluation  of  the  practicality  of  moment  methods 
using  standard  techniques.  The  long  run  times  for  a  mid¬ 


sized  scatterer  demonstrate  the  impractical ity  of  moment 
methods  for  dielectrics  using  standard  techniques. 


SCATl  ^.ING  OF  CYLINDRICAL  ELECTRIC 


FIELD  WAVES  FROM  AN  ELLIPTICAL 
DIELECTRIC  CYLINDRICAL  SHELL 


I.  Introduction 


Background 

The  radar  warning  receiver  (RWR)  on  the  F-16  is 
the  AN/ALR-69  built  by  Dalmo  Victor  Corporation.  Its  func 
tion  is  to  provide  warning  to  the  pilot  of  eneiry  radar 
activity.  It  informs  the  pilot  what  the  threat  is,  where 
it  is,  and  the  current  threat  status  (i.e.  search,  track, 
missile  launch,  etc.).  The  performance  of  the  ALR-69's 
direction  finding  (DF)  routine  is  degraded  by  a  side  or 
secondary  lobe  (SL)  located  30°  off  the  forward  position 
opposite  to  the  main  or  desired  lobe  [1] .  This  is  associ¬ 
ated  with  the  two  forward  antennas  only.  Figure  la  is  a 
sketch  of  the  desired  pattern.  Compare  this  to  Figure  lb 
which  is  a  sketch  of  the  actual  pattern.  Figure  2  shows 
the  location  of  the  antennas  on  the  aircraft.  The  SL 
causes  the  DF  routine  of  the  RWR's  processor  to  give 
erroneous  indications.  As  the  aircraft  approaches  the 
threat  emitter,  the  RWR  coit^ares  the  received  relative  sig¬ 
nal  strength  from  each  of  the  antennas.  The  threat  is  dis 
played  on  the  side  of  the  aircraft  that  received  the 


Left  Antenna 

Left  Antenna 

Desired  Antenna  Pattern 
(a) 

Actual  Antenna  Pattern 
(b)  [1:2] 

Figure  1.  Comparison  Between  Ideal  and  Actual 
Antenna  Patterns 
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stronger  signal.  However ,  due  to  the  secondary  lobe 
several  events  may  occur: 

1.  Both  antennas  receive  the  signal  equally  well; 
put  the  threat  in  front  of  the  aircraft. 

2.  Both  antennas  receive  the  signal  equally  well; 
the  RWR  can't  decide  where  to  put  the  threat  and  the 
result  is: 

a.  Displays  the  threat  location  correctly, 

b.  Displays  the  threat  on  the  wrong  side,  or 

c.  Displays  "flip-flops"  around  to  both  sides 
of  the  aircraft. 

3.  One  cuitenna  receives  the  signal  better;  dis¬ 
plays  the  threat  on  that  side  of  the  aircraft  which  may  or 
may  not  be  the  correct  side. 

Flight  tests  run  by  the  Tactical  Air  Warfare  Center  at 
Eglin  Air  Force  Base  have  verified  condition  3  in  that  the 
threat  was  displayed  on  the  wrong  side,  the  display 
"wandered"  around  on  the  screen,  and  then  jumped  to  the 
correct  side  as  the  aircraft  passed  by  the  target  [2] . 

Rationale 

This  situation  must  be  corrected  since  it  affects 
the  ability  of  the  F-16  to  perform  its  mission.  All 
threats  identified  forward  of  the  aircraft  are  suspect  as 
far  as  their  location  since  it  is  impossible  to  tell  if 
the  signal  is  being  picked  up  by  the  main  or  secondary 
lobe.  The  pilot  is  not  able  to  determine  vdiere  the  threat 


is,  and  he  has  no  idea  vdiat  maneuver  to  taUce  when 
approached  by  a  threat.  With  the  speed  of  the  F-16  and 
the  speed  of  an  approaching  aircraft  or  surface-to-air 
missile,  indecision  could  result  in  lost  aircraft  and 
lives.  It  should  be  noted  that  this  condition  has  been 
found  with  other  aircraft  RWR  systems  as  well.  Figure  3 
is  a  sketch  of  the  antenna  pattern  for  the  forward  RWR 
antennas  of  the  B-52  AN/ALR-46  [3] .  Figure  4  shows  their 
position  on  the  aircraft.  The  ALR-46  uses  the  same 
antennas  as  the  ALR-69. 

It  is  also  well  to  note  that  the  ALR-69  is  used 
on  the  A-10  aircraft.  In  this  case  there  is  no  side  lobe 


interference  with  radar  or  data  links)  [6]  or  biological 
tissue  (affects  of  microwaves  on  humans)  [7].  Holt, 
Uzunoglu,  and  Evans  were  concerned  with  the  meteorological 
scattering.  They  developed  eui  integral  equation  of  the 
Fredholm  type  directly  from  Maucwell's  equation  without 
reference  to  the  scatterer  [6] .  The  equation  developed 
was  a  general  second  kind  integral  equation  with  a  singu¬ 
lar  kernel.  The  singularity  was  removed  by  the  use  of  a 
transform  function.  The  result  is  a  pair  of  coupled  inte¬ 
gral  equations  that  were  solved  by  numerical  quadrature, 
vdiich  produced  a  numerically  stable  linear  algebra  equa¬ 
tion,  assuring  convergence. 

The  disadvantage  of  this  procedure  is  that  now 
there  are  two  equations  to  solve.  The  storage  require¬ 
ments  and  the  computational  time  increase  rapidly  as  the 
size  of  the  scatterer  becomes  large. 

Another  method  discussed  was  a  finite  element 
method  known  as  the  "Unimoment  Method"  [8;  9].  The  uni¬ 
moment  method  offers  the  eUaility  to  apply  the  radiation 
condition  without  the  use  of  complicated  programming  pro¬ 
cedures  as  required  in  the  past  use  of  finite  element 
methods.  Chang  and  Mel's  procedure  also  reduces  the 
storage  requirements,  when  compared  to  past  finite  element 
programs  [9:761] . 

Chang  and  Mei  essentially  solved  the  problem  using 
the  T-matrix  procedure.  After  creating  the  artificial 
boundary  outside  the  scatterer,  the  interior  problem  is 


6 


solved  using  the  finite  element  method.  The  field  between 
the  artificial  and  the  actual  boundary  is  approximated  by 
linear  combinations  of  functions  which  satisfy  the  Helm¬ 
holtz  wave  equation.  Chang  and  Mei  used  the  Fourier 
series  [8:36].  The  exterior  fields  were  expanded  in  terms 
of  Hankel  functions. 

While  this  procedure  does  offer  advantages  by 
separating  the  interior  and  exterior  problems,  for  large 
scatterers  a  significeuit  algebraic  equation  still  must  be 
solved.  Thus  this  method  will  only  handle  up  to  moder¬ 
ately  sized  scatterers  without  taking  considerable  time  and 
storage  resources. 

The  standard  method  used  in  the  past  has  been  the 
method  of  moments.  Richmond  [10]  used  this  technique 
in  solving  the  problem  of  scattering  from  infinitely  long 
dielectric  shells.  The  integral  equation  is  generated  by 
examining  the  polarization  current  that  results  when  the 
scatterer  is  illuminated  by  the  incident  field.  The 
unknown  electric  field  is  expanded  in  terms  of  pulse  func¬ 
tions. 

Since  Richmond  subdivided  the  shell  into  cells 
that  were  approximately  square,  integration  over  each  cell, 
including  the  self  cell,  could  be  accomplished  analytically 
The  cells  were  approximated  by  circles  with  a  radius  of 
half  the  length  of  the  cell. 

The  ability  to  analytically  solve  the  integration 
over  circular  cells  forces  the  shell  wall  to  be  thin.  If 


the  wall  is  thicker,  ass\iii^tions  will  have  to  be  made 
concerning  the  electric  field  in  the  shell.  Using  Rich- 
mond's  procedure,  the  electric  field  is  assumed  to  be  con¬ 
stant  over  a  cell.  For  anything  very  large,  this  would 
require  large  2unounts  of  storage  and  computational  time. 

It  is  possible  to  solve  this  scattering  problem 
in  terms  of  orthogonal  functions  such  as  Bessel  functions 
for  circular  cylinders  or  Legendre  polynomials  for 
spherical  objects.  In  the  elliptic  cylinder,  the  resultant 
functions  are  Mathieu  fiinctions.  There  is  not,  however, 
a  closed  form  solution  to  the  problem.  Dr.  Cavour  W.  H. 

Yeh  has  done  .considereible  work  with  scattering  and  travel¬ 
ling  wave  problems  in  connection  with  elliptic  shapes. 

In  each  case.  Dr.  Yeh  presented  a  solution  in  terms  of 
Mathieu  functions. 

In  a  study  of  sound  waves  scattering  from  pene¬ 
trable  objects,  pleme  waves  were  incident  on  to  an  ellipse 
at  different  aspect  angles  [11] .  The  resultant  patterns 
contained  side  lobes  that  would  move  as  the  angle  of  inci¬ 
dence  changed. 

The  primary  difficulty  with  this  type  of  approach 
is  the  Mathieu  functions.  Currently  a  library  does  not 
exist  on  the  base  cyber  computer  facility  for  generating 
Mathieu  functions. 


Previous  WorK 


Experimental  [1;  2],  General  Dynamics,  the  manu¬ 
facturer  of  the  F-16 ,  attempted  to  remove  the  side  lobe 
through  trial  and  error.  They  moved  the  antennas,  the 
radome,  added  material  to  the  radome,  the  airframe,  etc. 

It  was  through  this  work  that  some  important  information 
was  obtained. 

1.  The  SL  is  due  to  the  presence  of  the  radome. 
When  the  radome  was  removed,  the  side  lobe  disappeared 
[1:6]  . 

2.  The  SL  is  not  due  to  the  metal/dielectric 
boundary  between  the  airframe  and  the  radome.  The  radome 
was  moved  away  from  the  airframe  by  a  small  amount.  This 
introduced  a  new  boundary  layer  (metal  airframe,  air, 
radome  material)  and  the  lobe  became  larger  [2] 

3.  The  SL  is  not  being  diffracted  or  scattered 
by  the  dielectric/metal  boundary.  The  RWR  antennas  were 
moved  up  to  the  radome,  removing  the  break  between  the 
radome  and  the  airframe.  The  results  did  not  change  [2] . 
General  Dynamics  was  able  to  reduce  the  magnitude  of  the 
SL  by  adding  radar  absorbing  material  (RAM)  to  the  air¬ 
frame.  Figure  5  shows  the  amount  of  reduction  [1:5]  and 
the  location  of  the  RAM. 

Theoretical .  Schneider  [4]  modeled  the  radius  of 
curvature  of  the  radome.  He  used  a  circular  cylinder  large 
enough  so  that  the  curvature  of  the  cylinder  would  come 


close  to  the  curvature  of  the  F-16  radome  [4:5,35,37]  [5]. 
This  model  results  in  a  large  radius  (60  wavelengths) 
[4:2,35,37].  Schneider  used  two  methods,  a  series  solution 
and  a  numerical  approximation,  to  solve  the  boundary  value 
problem. 

In  either  method,  Schneider  was  xinable  to  repro¬ 
duce  the  SL.  He  was  hcunpered  by  a  small  computer  solving 
the  large  matrix  (3700  x  3700  elements)  that  resulted 
from  applying  point  matching  with  the  moment  method.  There 
fore,  he  presented  results  good  only  for  the  cylinder  with 
radii  of  0.6  and  6.0  wavelengths  [4:31,32;  Appendix  E] . 

The  series  solution  for  the  60  wavelength  scatterer  was 
calculated,  but  the  validity  of  the  results  is  question¬ 
able  since  deep  nulls  down  to  zero  appear  in  the  plots, 
and  the  nojnnalized  maximums  never  go  up  to  one  [4:29,30]. 
Schneider  admits  that  there  is  an  error  in  the  series 
solution  where  the  coefficient  of  one  term  is  half  of  the 
correct  value  [3:ii,iii].  This  error  will  effect  the  60 
wavelength  plot  and  may  account  for  the  deep  nulls  and  low 
peak  values. 

Proposed  Solution 

This  paper  will  present  a  discussion  on  the 
scattering  of  cylinderical  waves  from  a  large  dielectric 

^Note  that  this  differs  from  the  model  used  in 
this  thesis.  In  this  work,  the  radome  is  being  modeled 
as  an  elliptic  cylinder.  Schneider  modeled  the  curvature 
of  the  radome. 


scatterer.  The  scatterer  wxll  be  £ui  infinitely  long, 
cylindrical  shell  with  an  elliptic  cross  section  (see 
Figure  6).  This  is  2m  improvement  over  Schneider's  model 
in  that:  (1)  the  ellipse  represents  a  closer  approximation 
to  the  F-16  radome,  and  (2)  the  ellipse  has  a  much  higher 


a  -  semi-major  axis 

Py  b’ -  semi-minor  axis 

1  ^  c  -  focal  length 

Line!  Source 


o  o 


"o*S'^o 


Figure  6.  Dielectric  Scatterer 


rate  of  curvature  toward  the  tip  of  the  surface.  The 
increased  curvature  around  the  tip  of  the  ellipse  will 
cause  the  scattered  wave  to  be  "thrown"  in  a  particular 
direction.  This  hypothesis  was  suggested  by  the  work  done 
by  Chang  and  Mei  [8:41]  and  by  Poggio  and  Miller  [12:210]. 


To  solve  the  integral  equation,  the  method  of  moments  will 
be  used.  Harrington  [13],  and  Stutzman  and  Thiele  [14:306- 
372] ,  provide  discussion  for  application  of  the  procedure. 


Assumptions 

In  doing  this  work  it  was  assumed  that  the  radome 
could  be  modeled  apart  from  the  metal  airframe.  This 
assumption  was  suggested  by  the  General  Dynamics  study  [1] . 
It  was  assvimed  that  the  thickness  of  the  dielectric  shell 
is  thin  when  compared  to  wavelength.  This  will  simplify 
the  numerical  analysis  and  is  consistent  with  Schneider's 
work  [4:28]. 

Scope 

This  study  is  a  theoretical  study  which  is  almost 
entirely  removed  from  the  actual  problem  except  for  the 
dimensional  data  and  the  basis  for  the  models.  The  problem 
will  be  limited  to  the  two-dimensional  case  only.  This 
will  not  affect  the  results  since  the  actual  antennas  are 
coplanar.  There  will  be  no  attempt  made  to  solve  the  lobe 
problem.  The  goal  is  to  obtain  polar  plots  of  the  total 
electric  field  to  see  if  the  side  lobe  is  caused  by 
scattering.  The  real  need  is  to  understand  what  is  taking 
place  as  the  electromagnetic  fields  cross  the  boundary  of 
free  space  and  the  dielectric  shell. 


13 


II.  Development  of  a  One-dimensional  Fredholm 


Integral  Eguation  of  the  Second  Kind 

Solving  the  problem  using  moment  methods  requires 
the  derivation  of  an  integral  equation.  The  equation  used 
herein  was  developed  by  Richmond  [10:335-336],  extended  to 
elliptic  cylinder  coordinates.  Considerable  difficulty 
was  encountered  with  the  reduction  from  an  area  integral 
to  a  line  integral.  This  was  further  compounded  by  the 
singularity  of  the  kernel.  This  discussion  will  include  a 
description  of  the  model,  the  Richmond  integral  equation 
(IE)  in  elliptic  coordinates,  the  handling  of  the  singu¬ 
larity,  and  the  reduction  of  the  IE  to  a  line  integral. 
Appendix  A  provides  further  insight  into  the  coordinate 
system. 

Model 

Figure  6  is  a  diagram  of  the  elliptic  scatterer. 
The  shell  of  the  ellipse  is  0.05  wavelengths  thick.  This 
value  was  chosen  to  insure  applicability  of  the  thin  shell 
approximation  and  to  be  consistent  with  Schneider's  work 
[4:28].  In  elliptic  coordinates  the  thickness  must  be 
defined  as  a  dimensionless  quantity  to  be  consistent  with 
the  elliptical  coordinate  y.  Let  the  thickness  of  the 
shell  be  defined  as  t.  If  y^  defined  the  mid-radii,  then 
the  outer  wall  of  the  cylinder  is  y  +  t/2  and  the  inner 


is  -  t/2.  From  Appendix  A,  when  v  --  90®,  b  =  c  sinh  y. 
This  means  that  for  the  inner  and  outer  walls, 

b  -  T/2  =  c  sinh  (y^  -  t/2)  (2.1a) 

o 

b  +  T/2  *  c  sinh  (y^  +  t/2)  (2.1b) 

Solving  for  the  arguments  of  the  hyperbolic  sines  and  sub¬ 
tracting  (a)  from  (b),  the  result  is 

T  =  sinh"^  (b  +  T/2^_  sinh~^^  '  '^^-)  (2.2) 

c  c 

Expanding  the  arguments  of  (2.2)  in  a  Taylor  series  about 
T  and  taking  the  small  argument  approximation  of  the 
inverse  hyperbolic  sine. 


or 

T  -  T/a 
2  2  2 

since  c  »  a  -  b  and  a,  b,  and  c  are  large  when  compared 
to  T,  T,  and  X.  This  development  was  due  to  Johnson  [15] . 

When  making  reference  to  a  dimension  (i.e.,  a  or  b) 
it  will  be  done  with  reference  to  an  ellipse  drawn  through 
the  center  of  the  shell.  It  is  assumed  that  the  thickness 
is  constant  for  all  360®  and  that  the  focal  length  is  the 
same  for  the  inner,  middle  and  outer  radii.  The  error 


introduced  by  this  assumption  is  small  since  the  shell 
is  thin. 

To  insure  a  two-dimensional  case,  the  shell  and 
current  sources  are  infinitely  long  in  both  the  +/-z  direc' 
tions.  The  relative  permeability  is  constant  for  all 
positions  in  the  shell.  The  value  for  is  four.  This 
will  not  be  altered  in  the  insuing  programming.  The 
dielectric  is  perfect  and  no  losses  will  be  accounted  for. 
The  line  source  is  infinitesimally  thin  and  has  no  losses. 


Richmond  Integral  Equation 

The  equation  developed  by  Richmond  is  applicable 
since  the  scatterer  is  a  dielectric  cylindrical  shell. 
The  equation,  used  by  Schneider  for  moderately  sized 
scatter ers  [4:14,20],  is 


E^(x,y)  =  E(x,y)  + 


(e^-l)E(x',y') 


H^^^NkP)dx'dy'  (2.4) 
[10:336] 


where 

E  (x,y)  =  ^  +  E* 


r‘(x,y) 


-k  “l 
o 

4<jL)e 


H 


(2) 


(kp)  u. 


[16:224] 


The  superscript  i  and  s  signify  the  incident  and  scattered 
field,  respectively.  The  primed  terms  are  the  source 
terms  on  the  shell.  Due  to  the  shape  of  the  scatterer. 


the  use  of  the  more  widely  understood  polar  cylindrical 


coordinate  system  could  not  be  used.  It  was  too  cumber¬ 
some  to  describe  the  angle-dependent  radii  that  is  present 
in  an  ellipse.  Therefore,  the  Richmond  IE  was  converted 
to  the  elliptical  cylindrical  coordinate  system. 

Elliptic  Coordinates .  The  elliptic  cylindrical 
coordinate  system  describes  one  of  two  conic  sections, 
depending  on  which  of  the  two  coplanar  parameters  is  held 
constant.  As  Figure  7  indicates,  if  the  angular  coordin¬ 
ate,  V,  is  constant,  the  surface  is  a  hyperbolic  cylinder. 
If  y  is  constant,  an  elliptic  cylinder  results.  The  angle 
V  describes  the  angle  of  the  asymptote  of  the  hyperbola 
that  intersects  the  point  in  question.  It  can  easily  be 
shown  that  the  relation  between  the  polar  and  elliptical 
angles  is 

(j)  =  tan”^  [(b/a)  tan  v]  (2.5) 

The  range  of  v  is  from  0  to  360®  and  the  range  of 
y  is  from  0  to  ».  The  definition  of  y, 

y  =  tanh~^(b/a)  =  cosh~^  (a/c)  =  sinh”^  (b/c)  (2.6) 

blows  up  when  the  ratio  of  a  to  b  is  one  (or  c  is  0).  See 
Appendix  A  for  a  more  detailed  description  of  the  coordin¬ 
ate  system. 


The  differential  area  functions  for  elliptic  coordinates 
are 

^2  2 
.  cosh  y  -  cos  V  » 

dA  *  c^  •  j  2  2 

«  sinh^  y  +  sin^ 

Both  of  these  expressions  are  equivalent.  Using  (2.7)  with 
(2.4),  the  Richmond  IE  in  elliptic  coordinates  becomes 


2  2 


#(y,v)  =  E(y,v)  + 


(co^ 


jj  E(y',v')  Ocp) 

s' 

2 

y*  -  cos  V')  dy'dv* 


(2.8) 


2 

Note  that  the  c  term  from  the  differential  area  function 
has  been  pulled  out  of  the  integral  as  well  as  the  -  1 
term. 


The  Singular  Kernel 

Since  the  unlcnown  electric  field  is  both  part  of 
the  kernel  and  the  separate  function,  the  equation  is  a 
Fredholm  Integral  Equation  of  the  second  kind  (FIE  II) . 

The  equation  accounts  for  the  interaction  between  the  line 
source  and  a  point  on  the  shell,  euid  it  accounts  for  the 
interaction  between  that  point  and  the  rest  of  the  points 
on  the  shell.  There  are  essentially  two  coordinate  systems 
within  the  framework  of  the  problem.  Figure  8a  shows  the 
coordinate  system  for  which  the  line  source  is  considered 
the  source  point  and  the  shell  is  considered  the  observa¬ 
tion  point.  In  Figure  8b  the  point  indicated  in  Figure  8a 
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Line  Source 
(a) 


Shell  Interaction 
(b) 


- 1 _ I 

Figure  8.  Coordinate  Systems 

is  now  the  source  point  and  the  observation  point  is  any 
point  on  the  shell  including  the  source  point  itself.  As 
the  electric  field  from  the  line  sources  contacts  the 
scatterer,  a  polarisation  current  is  induced.  This  current 
reradiates  the  incident  field  as  the  scattered  field 
which  in  turn  reacts  with  the  other  points  on  the  ellipse. 
Each  point  on  the  scatterer  is  affected  by  two  electric 
fields;  the  incident  field  from  the  line  sources  emd  the 
scattered  field  from  the  rest  of  the  shell. 

The  problem  in  the  mathematical  emalysis  of  (2.8) 
is  the  consideration  of  the  interaction  between  a  point  on 
the  shell  with  itself.  This  point  on  the  shell,  known  as 
the  source  region  or  the  self  cell,  causes  the  Hankel 
function  of  (2.8)  to  blow  up  as  kp  goes  to  zero.  The  source 
region  interaction  with  itself  must  be  considered  for 


this  may  make  a  non-neglec table. contribution  to  the  total 
field. 


The  contribution  of  the  self  cell  is  done  by  the 
determination  of  a  correction  term.  This  correction  term 
accounts  for  the  electric  field  generated  at  the  source 
region.  Chen  discussed  the  source  region  in  terms  of  the 
tensor  Green's  function  [17:1201-1204] .  In  a  volume  the 
electric  field  is  determined  at  an  arbitrary  point  outside 
the  source  region  by 


E(r.)  =  I  G(r, ,r)  •  J(r)dv 


\diere 


G(rj^,r)  =  juy  [I  +  <|)(rj^#r) 

^o 


4>  (r,  ,r)  =  — e 
4ir  I  r£-r 


rjkJvrl 


G  is  the  tensor  Green's  function, 

I  is  the  unit  dyadic,  and 

4»  is  the  free  space  scalar  Green's  function. 
Accounting  for  the  source  region  (2.9)  becomes 


E(r  )  =  PV 
o 


G(r^,r) 


J(r)dv  +  E  (r  ) 
c  o 


(2.9) 


(2.10 


tdiere  the  "PV"  denotes  the  principal  value  (PV)  integral 
and  E  (r  )  is  the  correction  term  for  the  source  region. 


Figure  9.  Contour  Integration  About  B 

Principal  Value  Integral .  Before  discussing  the 
determination  of  the  electric  field  correction  term,  it 
would  be  prudent  to  define  the  principal  value  (or  the 
Cauchy  principal  value)  integral.  The  principal  value 
(PV)  integral  of  f(x)  over  [a,c]  when  f(x)  is  singular  at 
X  a  b  is  defined  as 

rc  lim  rb-r 

Pv/  f (x)dx  =  r -*•  ot  f (x)dx  +  I  f(x)dx]  (2.11) 

•Ja  Ja  Jbtr 

provided  the  limit  exists  [18:195] .  Essentially,  the  PV 
integral  is  eui  integral  taken  across  the  interval  excluding 
the  point  of  singularity.  At  that  point  a  contour  (or  area 
or  sphere)  is  taken  around  the  singularity.  The  radius  of 
that  contour,  r,  is  then  allowed  to  approach  zero  (see 
Figure  9) .  The  correction  term  is  the  field  from  across 
the  excluded  contour.  The  two-dimensional  case  involves 
a-  circular  area  and  the  three-dimensional  case  a  spherical 
volume  around  the  point. 


The  requirement  that  the  limit  exist  is  met  in 


the  case  of  the  Hankel  function.  According  to  A.  N. 
Tychonov,  the  improper  integral 


(2.12a) 


does  converge  as  long  as  «  < 3  [19:294].  In  the  two- 
dimensional  case,  the  integral  converges  for  «  <  2.  The 
Hankel  function  is  singular  as  In  |r|  .  However,  1/r  is 
more  singular  than  In  |r|.  Therefore,  the  area  integral 
of  the  Hankel  function  exists. 

Determination  of  the  Correction  Term.  The  value 
of  E^(r^)  is  a  function  of  the  shape  of  the  volume  excluded 

CO 

in  the  integral  evaluation.  In  calculating  E  ,  Chen  deter- 
mined  the  surface  charge  density  using  conservation  of 
charge.  The  electric  field  can  then  be  determined  by  the 
gradient  of  the  potential  due  to  the  surface  charge. 

E^(r^)  =  -V<t>  (2.12b) 


Yaghjian  presented  a  complete  discussion  of  the 
determination  of  the  Green's  function  in  NBS  Technical  Note 
1000,  A  Direct  Approach  to  the  Derivation  of  Electric 
Dyadic  Green ' s  Functions  [20] .  This  approach  does  not 
utilize  delta- function  techniques,  but  determined  a  general¬ 
ized  electric  dyadic  Green's  function  which  is  valid  in 
the  source  region.  Yaghjian  provided  a  table  of  the 


correction  terms  for  various  principal  value  geometries. 
He  determined  the  correction  factor  in  terms  of 


In  the  two-dimensional  case 


[20:56]  (2.13a) 

(See  Figure  10.) 

the  correction  term  becomes 

[20:59]  (2.13b) 

(See  Figure  11.) 


There  the  electric  field  is  now  determined  by  equation 
(2.14)  . 


E(r)  =  joju  lim  f  S  •  J  dv'  +  [20:12]  (2.14) 

^v-K)  Jv;-v  ^w^o 

e  J  e 

Or,  in  the  two-dimensional  case 


E(r)  =  W  r  G.JdA'+:^  [20:40]  (2.15) 

e  J  e 


In  (2.14)  and  (2.15),  the  V  -V^  and  A--A^  represent  the 
integral  excluding  the  source  region,  i.e.  the  principal 
value  integral. 

Tables  1  and  2  present  correction  terms  for  vari¬ 
ous  principal  geometries.  The  reader  is  encouraged  to 
read  NBS  Technical  Note  1000  for  further  insight  into  this 
problem. 


Table  1.  Tabulation  of  Source  Dyadic  L,  emd  Correspondence 

to  Previous  Authors  [20:60] 


PRINCIPLE  VOLUME 
(r  AT  CENTER) 

L 

PROBLEM  & 
AUTHOR(S) 

( 

SPHERE 

• 

s 

I 

T 

■Fm-?PJCE 

WILCOX  (1957) 
VANBLADEL  (1961) 

RIGHT  Cl 

RCULAR  C 

YLINDER 

®0 

(l-COS  60)6^ 

+  H  COS  0Q 

Cl 

A 

e. 

a, 

RCU 

i 

• 

a 

b- 

LAR  OR  S 
PENCIL 
f7| 

b 

►0  D- 

qUARE 

•  -»-  0 

-►0 

“t 

-r 

®x 

A 

✓ 

RECTANGULAR 

L  “ 

*  “yVy  *  “x'z'z’  <“x  *  «y  *  “z  - 
fiy,  and  0^  ARE  TWICE  THE  SOLID  ANGLE 

SUBTENDED  BY  A  SIDEXTO  THE  x,  y,  and  z, 

DIRECTION  RESPECTIVELY. 

i 

CUB 

E 

-► 

T 

3 

RECTANGULAR  CAVITY 

RAHMAT-SAMII  (1975) 

PILL  BOX 

h-H 

•>p 

6  6 

2  2 

RECTANGULAR  WAVEGUIDE 
AND  CAVITY 

TAI  (1973,  1976) 

WAVEGUIDE  ONLY 

COLLIN  (1973) 
RAHMAT-SAMII  (1975) 

Table  2.  Tabulation  of  2-D  Source  Dyadic  i  [20:61] 


In  the  two-dimensional  case,  if  the  electric  field 
is  aligned  with  the  vector  of  infinite  length  (i.e.  the  z 
axis  in  this  case) ,  then  the  electric  field  is  determined 
by 

-f  dA'  [21:261]  (2.16) 

and  the  correction  term  is  zero.  This  is  easily  seen  via 
Chen's  approach  since  the  surface  current  distribution  is 
zero  in  this  configuration  since  the  normal  is  orthogonal 
to  J.  For  example,  for  the  finite  cylinder,  the  current 
distribution  is  only  across  the  top  and  bottom  caps  of  the 
cylinder  [17:1203] .  In  the  two-dimensional  case,  the  top 
and  bottom  of  the  cylinder  are  at  +/-<»,  and  malce  no  con¬ 
tribution  to  the  field.  Loolcing  at  the  integral  that 
Richmond  evaluated  in  closed  form  for  the  self  cell 
[10:336];  for  the  case  where  n  =  m,  Richmond  determined 
that  the  self  cell  integral  of  the  Hanlel  function  was 


(ka)  adp'dij)'  =  j/2(TT]caH^^(]ca)  -  2j) 


(2.17) 


Richmond  had  approximated  the  square  cells  with 
circular  cells  of  radius  a.  If  a  is  allowed- to  go  to  zero, 
then  the  value  of  the  integral  goes  to  zero.  This  was 
verified  with  Bessel  tables  and  a  hand  calculator.  In 
Richmond's  formulation,  there  was  no  PV  integral; 
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therefore,  a  has  a  finite  value  and  the  source  region 
(m  =  n  cell)  had  to  be  calculated. 


It  would  be  well  to  note  that  not  all  of  the 
literature  agree  in  the  calculation  of  the  correction 
term.  J.  J.  H.  Wang  compiled  a  list  of  discrepancies 
[30:3-1] .  Table  3  and  Figure  12  list  his  results.  Since 
there  is  no  correction  term  involved  in  the  two-dimensional 
case,  this  does  not  impact  in  this  work. 


Reduction  of  the  Double 
Integral  to  One 

According  to  what  has  just  been  discussed,  the 
principal  value  integral  can  be  taken  over  the  shell  with¬ 
out  any  major  difficulty.  It  would  sin^lify  the  solving 
of  the  integral  equation  if  it  could  be  reduced . to  a  single 
integral . 

The  integral  in  (2.8)  is  repeated  in  expression 

(2.18)  . 


PV 


r  u, 

o  y^4/2 


+t/2  _ 

E(y’,v’)  H^'^'(kp) 


2  2 

(co^  y'  -  cos  v')  dy'dv' 


(2.18) 


where 


,  1  kV 

V-  4“ 


p  =  distance  function  between  two  points  in 
elliptic  coordinates 


See  Figure  12 


»  c  [  (cosh  u  cx>s  V  *■  cosh  y '  cos  v ' ) 

2  ^ 

+  (sijih  y  sin  v  -  sinh  y'  sin  v') 

Note  that  the  constant  "c"  is  the  focal  length  of  the 
ellipse,  not  the  speed  of  light.  The  "c"  was  chosen  since 
this  is  standard  notation  for  the  focal  length. 

Since  the  thickness  of  the  shell  wall  is  small 
(0.05  wavelengths)  it  was  assumed  that  the  electric  field 
was  constant  through  the  shell  thickness.  This  assumption 
is  backed  up  by  the  work  done  by  Lee,  Sheshadri,  Jamnejad, 
and  Mittra  [22] .  They  showed  that  through  dielectric 
shells,  with  a  Sickness  of  0.5  wavelengths,  that  the 
electric  field,  pattern,  was  only  slightly  altered  from  the 
no  dielectric  case  (see  Figure  13)  [22:3t7] .  Therefore 
with  T  being  X/20,  the  effect  of  the  shell  oh  the  field 
would  be  minimal.  The  electric  field  could  then  be  pulled 
out  of  the  first  integral. 

Figure  14  is  a  diagrcun  of  the  integral  around  the 
self  cell.  The  source  region  has  an  area  A--,  which  goes 

oK 

to  zero  as  y^  and  go  to  zero.  Since  the  integral 
exists,  the  limits  about  and  y^  may  be  taken  indepen¬ 
dently  of  each  other.  If  goes  to  zero  first,  the  areas 
indicated  by  a  "1"  become  very  small  and  their  total  con¬ 
tribution  may  be  neglected.  Therefore,  the  singularity 
becomes  a  line  singularity  and  the  principal  value  is 
taken  over  v  only. 


Figure  13.  E-plane  Radiation  Pattern  Through  a 
Spherical  Radome  [22:277] 


Since  the  electric  field  is  assumed  to  be  constant 


across  the  thickness  of  the  shell,  the  y  coordinate  is  a 
constant  and  may  be  integrated  directly.  This  procedure 
is  the  same  used  by  Stutzman  and  Thiele  in  the  derivation 
of  Pocklington' s  integral  equation  for  wire  scatterers 
[14:307-310].  Therefore,  the  resultant  integrand  is 

-cos^v')dv'  (2.19) 
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H^(kcPo) 


CO^  u 


The  Richmond  integral  equation  reduced  to  one  variable  is: 


_  — 

?  (y,v)  =  E(u,v)  +  Xp  tPV  J  E(v1q,v 


•) 


•  H  ^^^(kcp J  (co^^  u  -  cos^  V')  dv'  (12.20) 


where  (from  2.4) 


-k  ^  I 

^  •H''"  -  Ho  ’  “z 

o 


Xp  -  j3k2c2/4 


p  »  V(oo*uoos^-ooahu  cos\; +(sinhpsii^;-sinhii  sinv')^ 


In  doing  this  integration,  it  was  assumed  that  the 
focal  length,  c,  is  constant  for  the  entire  shell.  This 
is  justified  by  the  thin  shell  approximation.  Therefore, 


by  the  Method  of  Moments 

Method  of  Moments 

General  discussion  on  moment  methods  and  applica¬ 
tion  to  wire  antenna  problems  may  be  found  in  any  of  the 
newer  antenna  texts.  Stutzman  and  Thiele  [14:Chapter  7] 
provide  2m  excellent  tutorial  for  the  new  student.  Unfor¬ 
tunately#  discussion  on  the  use  of  moment  methods  for  non- 
metallic  scatterers  is  severely  limited.  Mittra's  book. 
Computer  Techniques  for  Electromagnetics,  has  a  discus¬ 
sion  of  such  problems  in  a  very  general  sense  [28]. 
Harrington  [13]  discusses  the  dielectric  problem;  however, 
the  discussion  is  based  exclusively  on  Richmond's  results 
[10]  . 

Moment  Methods  are  a  general  procedure  converting 
problems  posed  as  integral  equations  into  linear  equations, 
which  then  can  be  solved  numerically.  The  method  of 
moments  is  a  projection  method  in  which  the  unknown  is 
approximated  or  projected.  The  approximated  equation  is 
then  solved  exactly  [23:3].  Consider  the  relation 

L(f)  -  g  (3.1) 

where  L  is  a  linear  operator  operating  on  f,  and  g  is  a 
known  result.  The  unknown  to  be  determined  is  f.  Expand 


£  by  a  series  of  functions  in  the  domain  of  L  such  that 
N 

f»Zoc  f  (3.2) 

n  " 

These  functions  are  known  as  basis  or  expansion  functions. 
Letting  L  operate  on  £,  (3.1)  becomes 

N 

Z«nL(f^)=g  (3.3) 


An  inner  product  is  then  tcUcen  with  a  weighting  or 
testing  function  and  L(f^)  and  g,  respectively.  This  func 
tion  is  in  tlie  range  of  L.  Therefore  3.3  becomes 

N 

Z  «  <w  /L(f  )>  *  <w  .g>  m=  1,2,3,...  (3.4) 

^  n  m  n  m 

n 

In  matrix  form,  this  is 


'U  '-n'  -  <V 


where 


im 


<Wj^,L(fj^)> 

<W2,L(fj^)> 


<W2,L(fj^)> 


<Wj^,L(f2)> 

<W2rL(f2)> 

<W2,L(f2)> 


<Wj^,L(f2)> 

<W2,L(f2)> 

<W2,L(f2)> 


(3.5a) 


<  Wn'I'iy 


Solving  for  with  a  nonsingular  matrix,  the  unknown  f  can 

be  determined.  The  selection  of  f  and  w  is  as  much  art 

n  m 

as  science  and  has  been  widely  discussed  in  the  litera¬ 
ture.  A  thorough  treatment  will  be  given  later  in  this 
chapter.  This  discussion  on  the  general  theory  of  moment 
methods  is  from  Harrington  [13:3-7]. 

With  this  brief  introduction,  the  intent  of  this 
chapter  is  to  apply  the  method  of  moments  to  solve  a 
Fredholm  Integral  Equation  of  the  second  kind;  specifically 
(2.20).  This  chapter  will  include  the  theory  of  moment 
methods  and  the  application  to  the  elliptic  scattering 
problem. 


General  Application 

Starting  with  equation  (2.20), 


E^(y,v)  «  E(y,v)  +  E(yjj,v' 


) 


•  (kcP)  (cosh^  y„  -  cos^  V)  dv'  (3.6) 


(2) 


expand  the  unkno\m  electric  field  on  the  dielectric  shell 


in  terms  of  the  general  basis  function 

E,  -  2  C  P(\> -V)  (3.7) 

z  ~  n  n 

(Note  that  in  (3.6)  the  "cut  integral"  symbol^,  is  used 
instead  of  "PV" .  This  signifies  that  the  singularity  is 
not  included  or  has  been  "cut  out").  This  makes  (3.6) 
become 

i 

E2(y,v)  =  Z  C^P(Vj^-v)  +  V  T  ^ 

Jo 

*  H  ^(kcp)  (cosh^  y  “  cos^  v')  d\;'  (3.8) 
o  o 

Now  apply  the  testing  function  to  (3.6)  and  (3.8). 
Moving  the  summation  sign  outside  of  the  integrals  and 
taking  the  inner  product  (as  shown  in  (3.4),  we  obtain 


r^nH-1 

fU^*r/2 

E^(yrV,y  ,v  )  W(v  -v)  (cosh^  p  -  cos^  v)  dpA; 
s  s  m 

2 

c 

J 

'^m-l 

J 

U^-T/2 

i^C 

n 

r^mhl 

J 

^m-1 

J 

U„-r/2 

W(v  -v)  P(v  -v)  (cosh^p  -  cos^  v)  dpA;  (3.9 
m  n 

r 

.Vl  ^ 

*  “n  j 

J  t 

Vl  y  -r/2  ° 

(2)  2  2  2  2 
*  Hq  (kcp)  (co^  p  -  cos  v)  •  (oodi  p^-cos  v')dv'dpdv 

it^*l  ,2,3. . . 


The  order  of  integration  may  be  interchanged  except  that 
the  integration  over  v '  must  be  done  before  the  integra¬ 
tion  over  V .  This  is  due  to  the  fact  that  v '  is  over  the 
entire  shell  and  is  part  of  the  IE  to  be  solved,  whereas 
the  integration  in  v  is  done  only  over  the  domain  of  the 
weighting  function.  This  is  not  the  most  general  since 

it  is  assumed  that  P(v_)  and  W(v  )  are  sub-domain  functions 

n  m 

Therefore  the  functions  are  zero  for  v  not  within  the 
domain  of  P(v) .  The  use  of  entire-domain  functions,  which 
are  valid  across  the  entire  range  of  v  (the  entire  scat- 
terer)  is  not  considered  here. 

In  this  discussion,  it  should  be  noted  that  a  two- 
dimensional  application  of  the  method  of  moments  is  pos¬ 
sible.  Essentially,  in  assuming  that  the  y  coordinate  was 
a  constant  and  could  be  integrated  out,  y  was  expanded  by 
a  pulse  basis  function  valid  across  the  entire  domain  of  y. 
The  integration  over  y  was  done  in  the  seune  manner  for  the 

voltage  term  (the  term  containing  E^) .  The  "middle" 

z 

integrals  (i.e.,  they  are  physically  located  in  the 
middle  of  (3.9))  were  done  directly  since  these  integrals, 
involving  only  the  differential  area  function,  exist  in 
closed  form.  After  the  integration  over  y,  (3.9) 


becomes 


3 


N 


2 

C  T 


"nrt-l 


f  2  2 

r|  E  (yQ»v,y„,v^)  W(v^-y)  (cosh  y^  -  cos  v) 


s  s  m 


'm-l 


I,  r 

^  [cosh  (2y  J  sirh(T)  +  t]  E  c^  I 
£  o  n  I 


Wi 


W(v  -v)P(v  -V)  dv 
in  n 


m-1 


c^tEc 


.  / 


Wl 


W(v  -v)P(v  -v)  cos  V  dv 
m  n 


(3 


Vl 


2  2 

+  A^tVec 


Vl  ,Vi 


W(v  -v)P(v  -V') 
m  m 


V  ,  "V  . 

m-1  n-1 


<21  2  2  2  2 
*  ’  (kp)  (codn  y^  -  cos  v’)  (com  y^  -  cos  v)  dv’dv 


where 


_k 

E^(H,,v,Ug,Vj)  -  (kp^)  \ 

o 

Xp  =  jSk^c^/A 

2  2 
p  *  [cosh  y^(c  cos  v  -  c  cos  v  ) 
s  o  s  s 

2  2  k 

+  sinh  y  (c  sin  v  -  c  sin  v  )  ]^ 
o  s  s 


y  ,v  ^  line  source  coordinates 
s  s 


y_,v,v'  ~  scatterer  coordinates  (primed  denotes 
°  source  region) 


.10) 
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Equation  (3.10)  can  be  written  as 


V  =  C  [Z  +  Z* 
m  n  m 


nn 


=  C  Z 


n  nn 


(3.11) 


where  the  vector  in  the  brackets  Z_,  is  the  single  inte¬ 
rn 

gral  from  the  total  electric  field  and  Z'  ^  is  the  double 

mn 

integral  over  v  and  v'  involving  the  Hankel  function. 
Together  these  make  up  the  reaction  matrix,  Z 

mn 

Solving  (3.11)  for  and  using  (3.7),  the  electric  field 
and  the  induced  current  can  be  calculated.  The  radiation 
pattern  can  then  be  determined  in  the  far  field  using 

Az(p)  = 

Ez  =  -jwy^Az  [14:25]  (3.12b) 


N 

E 


e'^oP 


itFl  VWv  ceU 


ff  ♦''as. 


m 


[16:229] 


(3.12a) 


where  based  on  the  integration  over  the  cell  of 

the  basis  functions  and  the  Hankel  function.  Equation 
(3.12c)  i.:  derived  with  H  =  V  x  A  consistent  with 
Harrington  [16:77]. 


Basis  Functions 

The  rest  of  this  discussion  will  center  on  the 
basis  and  weight  functions,  how  they  are  applied  in  moment 

^The  summation  over  n  from  1  to  N  shows  that  the 
cor  ribution  from  each  cell  must  be  summed  together  to 
obtain  the  total  vector  potential  at  a  point. 


methods/  and  the  set  of  functions  used  in  this  problem. 
Harrington  states  that  one  of  the  main  tasks  in  using 


this  problem  solving  technique  is  the  choice  of 

and  W(v  )  (P(v  )  and  W(v  )  are  the  basis  and  testing 
m  n  m 

functions,  respectively)  .  needs  to  be  as  close  an 

approximation  to  f  as  possible  and  linearly  independent. 
W(v^)  should  also  be  linearly  independent  and  chosen  so 
that  the  inner  products  with  g  depend  on  the  properties  of 
g  (the  solution  of  L(f)).  Additional  factors  are  the  solu¬ 
tion  accuracy,  relative  ease  in  the  evaluation  of  the 
matrix  elements,  the  number  of  segments  required  for  con¬ 
vergence,  and  the  stability  of  the  resultant  matrix  equa¬ 
tion  [13:7] . 

Below  is  a  list  of  the  standard  basis  and  weight 
functions  commonly  used.  These  are  sub-domain  functions 
since  they  are  defined  to  be  zero  outside  the  domain  of  the 
function. 


Piecewise  Uniform  (Pulse  Function) : 


Jj(Z) 


I. 

D 


<  Z  < 


0 


elsewhere 


(3.13a) 


Triangle  Fvinction  (Piecewise  Linear) ; 


lAZ  -  ZJ 


Z.  -  Z.  , 
3  3-1 


Z  ,  <  Z  <  Z. 
3-1  -  -  3 


I.(Z.  1  -  Z) 


Piecewise  Sinusoidal: 


(3.13b) 


Jj(Z)  = 


sin  [k(Z-Z^_^)] 
six  [)c(Z^  -  Zj_3^)] 

sin  (k(Z^^j^-Z)] 
sin  [k(Z^j^  -  Z^)] 


z .  .  <  z  ^  z . 
3-1  -  .  -  3 


Z.  <  Z  <  Z.^, 
3  -  -  3+1 


(3.13c) 


Suadratic  Interpolation: 


elsewhere 


Jj(Z) 


A.+  B.(Z-Z.)  +C.(Z-Z.)‘ 
3  3  3  3  3 


Zj  <  2  < 


elsewhere 


(3.13d) 


Sinusoidal  Interpolation : 


Jj(Z) 


Aj  +  Bj  sin  (k(Z-Zj)] 
+  Cj  cos  [k(2-Zj)] 


Z .  <  Z  <  Z 
3  -  -  3+1 

(3.13e) 

elsewhere  [4  7:23,24] 
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Figure  15  illustrates  three  of  these  functions. 

The  use  of  the  piecewise  functions,  forces  the 
partitioning  of  the  segments  into  two  regions  ranging  from 
to  Zj  and  from  Z^  to  The  result  of  the  integra¬ 

tion  is  summed  together  to  obtain  the  total  contribution 
for  each  cell.  This  division  is  due  to  the  derivative 

discontinuity  at  Z  [24:535].  Butler  and  Wilton  discuss 

n 

methods  of  modifying  certain  functions  in  order  to  increase 
the  convergence  rate.  Specifically  with  (3.13e)  in  one 
procedure  the  constants  and  are  adjusted  to  force 
the  current  and  its  derivative  to  be  continuous  at  each 
of  the  cell's  endpoints.  This  results  in  a  basis  set  vdiich 
causes  the  current  and  its  derivative  to  be  continuous 
across  the  cell.  In  another  procedure,  the  current  in  the 
jth  cell  is  required  to  satisfy  Equation  (3.14).  This  is 
known  as  extrapolated  continuity  [24:535]. 


(3.14) 


45 


Richmond  used  square  pulses  in  calculating  the 
field  in  the  dielectric  shell  [10:336] .  Use  of  such  func¬ 
tions  assumes  a  constant  field  across  the  shell,  which  is 
a  reasonable  approximation  if  the  shell  is  thin  and  the 
scatterer  small  (small  compared  to  wavelength) .  Hagmann, 
et  al.  showed  that  larger  cells  may  be  used  if  variations 
are  allowed  for  the  field.  In  the  plane-wave  correction 
method,  the  field  is  represented  by  a  supe imposition  of 
plane  waves  (3.15). 

E„(r',0')  =  Z  B.e~^’  ^  [25:744]  (3.15) 

z 

In  the  cy Under ical-cell  correction  method,  the  square 
cells  are  replaced  by  circular  cells  as  Richmond  did 
[10:336,337].  The  circular  cylinder  with  TM  excitation 
will  have  fields  with  the  cell  determined  by 

00 

E,  =  Z  b  J  (kr)  cos  (ne^C  )  (3.16) 

z  n  n  n 

where  b  and  C  are  determined  by  the  incident  wave 
n  n 

[25:746]. 

In  another  example  of  methods  for  improved  conver¬ 
gence,  Blue  used  different  basis  functions  depending  on 
the  geometry  of  the  problem  and  the  location  of  the  seg¬ 
ment  in  relation  to  that  geometry.  Blue  showed  that  a 
conducting  strip,  6 OX  wide,  could  be  done  with  only  17 
basis  functions  [26:1894].  Analyzing  the  problem  before 


handf  Blue  used  three  basis  functions,  each  determined  by 
the  domain  of  the  function  [26:1902].  As  an  example  with 
the  problem  being  discussed  in  this  paper;  for  the  more 
pointed  end  of  the  ellipse  (as  y  goes  to  zero) ,  Hankel 
functions  could  have  been  used.  This  is  based  on  the  func¬ 
tion  which  describes  the  electric  field  distribution  in  a 
curved  optical  waveguide  [27:2125,2126].  A  sinusoidal 
function  could  have  been  used  in  the  smooth  part  of  the 
ellipse  (as  x  goes  to  zero)  since  the  structure  approaches 
a  dielectric  slab  waveguide.  This  would  mean,  however,  a 
tradeoff  is  being  made  for  accuracy  and  convergence  rate 
versus  simplicity  of  implementation. 

Testing  Functions  and 
Solution  Methods 

Testing  Solutions.  The  choice  of  testing  functions 
is  the  same  as  for  the  basis  functions.  As  discussed 
earlier,  the  ideal  testing  functions  would  result  in  easy 
in^lementation ,  high  speed  computation,  and  a  fast  con¬ 
vergence  rate  to  produce  a  highly  accurate  solution. 
However,  it  is  not  possible  to  have  all  four  criteria  and 
be  able  to  solve  all  of  the  problems  encountered. 

Solution  Methods .  Table  4  provides  five  solution 
methods  or  function  pairs  commonly  used.  Each  of  the 
methods  of  Table  4  has  advantages  and  disadvantages.  Point 
Matching  provides  a  highly  accurate  solution,  but  the 
number  of  cells  required  maJces  this  functional  pair 
impractical  for  scatterers  of  any  significant  size.  The 


Table  4 .  Moment  Method  Functional  Pairs 


nth  Term  of  A 

Method 

(basis) 

Testing  Function 

Galerkin  a_b_(x) 


Least  Square  (x) 


Point  Matching  a^  <S  (x-x^) 


Point 

Collocation 


Subsectional 

Collocation 


a  b  (x) 
n  n 


6 (x-x  ) 


6 (x-x  ) 
m 


V(x  )  Z  a_  b  (x)  6 (x-x  ) 

n  -  no  D  m 


Q(x)  IS  a  positive  definite  function  of  position 
6(x)  is  the  delta  function  [12:188] 


chief  advantage  of  point  matching  is  the  easy  implementa¬ 
tion  and  quick  calculation  of  matrix  elements  (low  CPU 
time) .  Galerkin  functional  pairs  converge  faster  than 
least  square  pairs ,  but  a  Galerkin  pair  will  not  always 
converge  to  a  solution.  Sarkar  showed  that  the  least 
squares  method  will  always  produce  meaningful  solutions , 
even  when  the  solution,  g,  is  not  within  the  range  of  L 
(the  linear  operator  from  (3.1))  [23:2].  The  results  of 
this  study  showed  that  Galerkin 's  method  also  takes  con¬ 
siderable  CPU  time  to  calculate  matrix  elements  when  the 
elements  have  to  be  integrated  numerically.  This  statement 
would  be  true  for  any  of  the  methods  involving  functions 
that  could  not  be  integrated  analytically. 


9 


After  careful  consideration,  the  method  used  to 


solve  the  dielectric  elliptical  shell  problem  was 
Galerkin's  method  with  the  piecewise  sinusoidal  function 
set.  This  functional  set,  of  the  established  ones,  most 
closely  resembled  the  electric  field  as  it  propagates  in 
a  dielectric  slab.  The  use  of  Hankel  functions,  while 
suggested  by  Blue's  article,  would  be  too  conqplicated  to 
use  even  if  a  better  match  is  obtained.  Other  reasons  for 
choosing  this  set  were:  one  of  the  most  popular  methods, 
well  documented,  and  ease  in  understanding  implementation. 
The  biggest  drawback  is  the  need  for  numerical  integra¬ 
tion  for  each  matrix  element. 


Therefore,  the  basis  and  testing  functions  are: 


(Sin  +,-v)]  sin  [k(v'-v  ,)] 
E(yQv')  ~  'Cj^  jsin  sin 


(3.17) 


Please  note  that  the  (v-v')  of  equation  3.17 
is  the  distance  between  the  angles  at  v  and  v ' .  The  dimen 
sion  is  in  meters.  This  is  consistent  with  the  fact  that 
the  dimention  of  the  wavenumber,  k,  is  rad/m. 

A  good  reference  for  the  implementation  of  the 
simusoidal  basis  function  is  Richmond's  report.  Computer 
Analysis  of  Three-Dimensional  Wire  Antennas  from  The  Ohio 


State  University  [29]  .  While  dealing  with  metallic 
scatterers,  the  features  such  as  segment  division,  inte¬ 
gration,  etc.,  are  demonstrated  better  than  in  the  text 
book  dipole  implementation. 


I 


k-' 


IV.  Programming  the  Moment  Method  Solution 


The  final  step  in  the  procedure  is  the  generation 
of  code  to  implement  the  method  of  moments  and  solve  the 
linear  algebra  equation,  (A)x  =  b.  The  Fredholm  integral 
equation  of  the  second  kind. 


E^(y,v)  =  E(u,v)  +  \pT  -L  E(u^,v)  (kcp) 


2  2 

•  (cosh  -  cos  v')  dv* 


(4.1) 


has  been  reduced  to  a  linear  algebra  problem.  The  matrix, 
A,  consists  of  two  elements,  and  where 


ymH  ^'^mtl 


/  P(v„-V)  W(v^-v) 


Vi  Vi 


.  H  (kcp)  (codi^  V  -  cos^v') 
o  o 

2  2 

•  (cosh  u  -  cos  v)  dv'dv 
o 


(4.2a) 


,  ,Vl 

/  / 

vi  V/2 


P(v  -v)  W(v  -v) 
n  m 


2  2 

(cosh  y  -  cos  v)  dydv 


(4.2b) 


2Uid  the  voltage  vector  is 


''m- 


.22 

■*0  = 

4E»e 


/ 


mH 


P(v„-v)  H 
m  o 


(2) 


o  s 


itt-1 


(4.3) 


•  (cosh^  -  cos^) 


where  k  is  the  wave  number  in  the  dielectric  if  integrating 

for  or  /y^e^  for  free  space  when  calculating 

Vmxt  P(^  and  W(v  -V)  are  the  sinusoidal  basis  func- 

nM  n  m 

tions  just  discussed. 

To  calculate  the  matrix  elements  and  solve  this 
the  linear  algebra  problem  will  require  special  function 
calculators,  numeric  integrators,  and  linear  algebra 
solvers.  The  self  cell  (cell  where  n  =  m)  will  require 
special  handling  as  will  the  cells  at  the  +/-90°  points. 

This  chapter  will  discuss  the  programs  developed  to  solve 
the  linear  algebra  equation.  The  discussion  will  include 
a  description  of  the  special  functions  euid  routines.  The 
chapter  will  close  with  a  discussion  on  the  problems  still 
unresolved  with  the  program. 

Machine 

Due  to  the  large  core  and  time  requirements,  this 
program  was  processed  on  the  Aeronautical  Systems  Division's 
CYBER  computer  system.  The  ASD  system  consists  of  two 
Control  Data  Corporation  (CDC)  mainframes,  CYBER  74  and  750, 
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which  operate  In  parallel.  The  750  has  262,000  words  of 
central  memory.  The  CYBER  74  has  131,000  words.  The 
CYBERs  have  14  significant  digits  for  real  variables.  The 
machines  support  a  variety  of  languages  and  support 
packages.  The  system  also  includes  interactive  processing 
(INTERCOM) ,  plotters  (CALCOMP  and  DISSPLA) ,  and  special 
libraries  (IMSL,  FUNPACK,  etc.)  [31:3]. 

Proqrzun  Overview 

The  language  used  in  this  program  was  FORTRAN  V, 
with  COC  extensions.  Version  5  complies  with  the  American 
National  St2Lndards  Institute  FORTRAN  77  [32:V]  .  The  CDC 
extensions  used  were  minor  consisting  of  the  use  of  sine, 
cosine,  and  tangent  functions  in  degrees,  hyperbolic 
arctangent,  the  CPU  SECOND  functions  (returns  current  pro¬ 
cessing  time  elapsed  since  start) ,  and  CDC  FORTRAN  control 
statements.  Due  to  the  heavy  use  of  the  IF-THEN-ELSE 
statements  and  the  use  of  zero  index  values,  this  progreun 
will  not  compile  on  a  FORTRAN  4  compiler  without  major 
modifications.  The  programs  were  nm  with  the  optimizer 
set  to  1  (0PT=1) ,  vdiere  the  compiler  will  optimize  the 
code  by  the  following  steps: 

1.  Redundant  instructions  and  expressions  are 

removed . 

2.  PERT  critical  path  scheduling  is  done  to 
utilize  the  multiple  functional  units  efficiently. 


3.  Svibscript  calculations  are  sinplified  euid 
values  of  simple  integer  variables  are  stored  in  machine 
registers  throughout  loop  execution  for  certain  do  loops 
[32:11-7,8]. 

Despite  the  reduction  in  the  amount  of  cells 
required,  the  storage  requirements  for  the  progreun  were 
large  when  the  full  sized  scatterer  (semi-major  cucis  of 
12X,  semi -minor  axis  of  6X)  was  run.  It  was  not  possible 
to  run  all  of  the  algorithms  in  one  job.  Therefore,  it 
was  necessary  to  break  up  the  program  into  three  programs 
which  were  completely  separate.  Figure  16  has  a  flow 
chart  of  the  simple  structure.  The  TRANSF  control  card 
was  used  to  control  the  flow.  Using  the  TRANSF  card,  the 
input  program  ran  first,  then  the  next  job  listed  on  the 
TRANSF  card  ran.  Program  3  could  not  run  until  2  was 
finished  and  2  could  not  run  until  the  first  program  was 
completed.  This  offered  several  advantages.  Each  module 
could  be  developed  separately.  Since  the  input  program 
created  a  perm^ulent  file,  once  the  geometry  was  set  up  and 
the  file  created,  it  did  not  have  to  be  mn  again.  It  was 
possible  to  reduce  core  requirements  which  decreased  turn¬ 
around  time  by  running  each  program  separately.  This 
structure  was  much  easier  to  understand  and  develop  than 
using  segmentation  or  overlays. 
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Figure  16.  Job  Stream  on  CYBER  175 
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PROGRAM  3 
CALCULATES  TOTAL 
FAR  FIELD.  PLOTS 
RESULTS  USING 
DISSPLA 


Figure  16 — Continued 


Library  Subroutines 

Due  to  the  nature  of  the  equation  being  solved, 
special  functions  and  procedures  were  needed.  These 
included  elliptic  integrals,  Bessel  functions,  integrators 
(numeric  quadrature),  and  linear  algebra  equation  solvers. 
The  ASD  computer  system  has  many  functions  and  routines 
already  built  in  through  several  different  program  packages 
Rather  than  go  through  the  procedure  of  developing  and 
testing  these  programs,  the  "ccmned"  programs  were  used. 
These  offer  the  advantage  that  they  are  well  known  through¬ 
out  the  industry  and  tried  and  tested  on  many  different 
machines. 

Functions.  The  complete  elliptic  integral  of  the 

second  kind  was  used  to.  calculate  the  circumference  of  the 

ellipse.  The  Bessel,  J  (x) ,  and  the  Neumann,  N_(x)  (or 

o  o 

Y  (x) ) ,  functions  were  used  to  calculate  the  matrix  ele- 
o 

ments.  The  library  selected  to  calculate  these  functions 
was  the  FUNPACK  library  [33] .  FUNPACK  is  a  library  of 
functions  and  subroutines  that  can  return  the  results  of 
Bessel  fxinctions  of  first  euid  second  kind,  modified  Bessel 
functions,  elliptic,  exponential,  Dawson's  integrals,  and 
other  functions.  There  are  13  programs  in  the  current 
library  which  was  released  in  1976  as  FUNPACK  Release  2. 
FUNPACK  was  developed  as  part  of  the  National  Activity  for 
Testing  Software  (NATS)  project.  Obtained  from  the  Argonne 
Code  Center,  FUNPACK  was  specifically  designed  for  the 
CDC  6000-7000  machines.  The  documentation  that  accompanied 


ti 


the  program  [33] ,  provides  the  description  on  how  to  use 
the  code,  the  accuracies  and  limitations  of  the  routines, 
and  where  they  were  tested.  The  results  of  the  Bessel 
functions  for  euid  were  heuid-checked  against  tables 
and  found  to  be  accurate.  The  only  limitation  to  accuracy 
is  the  host's  ability  to  accurate  the  values  of  the  basic 
functions  such  as  natural  log,  sine,  cosine,  etc.  [33:11]. 
The  nicest  part  of  the  FUNPACK  library  is  that  the  package 
was  developed  for  the  CDC  processors  and  take  full  adveui- 
tage  of  the  14  working  numbers  of  accuracy  in  standard 
precision.  It  was  for  this  reason  that  the  FUNPACK  library 
was  chosen  over  the  better  known  International  Mathematical 
and  Statistical  Library  (IMSL) .  It  should  be  noted  that 
the  FUNPACK  libraries  are  still  being  supported  by  Argonne 
and  will  continue  to  receive  support.  Changes  will  be 
supplied  automatically  to  users  [33:11,22-23].  Another 
adveuitage  of  the  library  is  that  the  function  calls 
resemble  the  name  of  the  function,  making  it  much  easier 
for  the  future  user  to  read  the  program  and  understand 
what  operation  is  being  performed  without  a  great  number  of 
comment  statements. 

Integrators.  As  can  be  seen  from  equations 
(4.1-3),  the  numerical  integration  is  an  important  part  of 
the  program.  The  integrals  of  (4.2a)  and  (4.3)  do  not 
exist  in  closed  form.  While  the  single  integral  of  (4.2b) 
can  be  done  analytically,  the  integral  over  the  angle,  v, 
was  done  numerically  to  reduce  the  chance  of  algebraic 
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error.  The  radiation  routine  used  to  calculate  the  far 


field  also  needs  a  single  integrator.  The  IMSL  has  a 
single  and  double  integrator.  The  documentation  on  IMSL 
is  limited  and  the  accuracy  of  the  routines  is  not  known 
except  to  rely  on  the  reputation  of  the  library.  The 
AFIT  Digital  Computer  Manual  called  the  package  "...  the 
recommended  routines"  [34:55]. 

The  single  integration  routine,  DCAORE,  uses  a 
cautious  Adaptive  Romberg  Extrapolation  Algorithm.  DCADRE 
is  computed  as  the  stun  of  estimates  for  the  integral  of  a 
function,  F(x),  over  suitably  chosen  subintervals  of  the 
limits  of  integration.  If  the  routine  is  unable  to  find  an 
acceptable  estimate  on  a  given  subinterval,  the  subinterval 
is  divided  and  each  of  the  new  subintervals  is  handled 
separately  [35:DCA0RE>1, 2] .  It  is  because  of  this  process 
of  interval  subdivision  that  DCADRE  can  take  a  great  deal 
of  time.  The  acceptability  of  an  interval  is  determined 
by  a  relative  error  input  in  the  function  calling  state¬ 
ment.  If  the  routine  can  not  get  within  the  limits  sup¬ 
plied,  then  the  function  writes  out  an  error  message  and 
supplies  the  best  answer.  However,  even  if  the  error  mes¬ 
sage  is  written,  the  best  answer  is  still  better  than  the 
standard  integration  routines.  DCADRE  may  return  wrong 
answers  if  the  frequency  of  the  integrand  is  very  high, 
but  this  problem  may  be  overcome  by  dividing  the  interval 
and  calling  the  routine  several  times  [35:DCADRE-2] . 


The  double  integrator  is  DBLIN,  which  uses  DCADRE 
to  do  the  single  integration.  However,  use  o£  the  Romberg 
integration  routine  in  two  variables  increased  run  times 
significantly.  Therefore,  the  integral  over  was  done 
via  a  Sin^son  rule  integrator  [38:311]  while  DCADRE  inte¬ 
grated  over  V  . 

n 

Matrix  Equation  Solvers .  The  IMSL  has  two  routines 
for  the  solving  of  complex  linear  algebra  problems:  LEQTIC 
and  LEQ2C.  There  is  also  a  library  of  linear  algebra 
routines  from  the  Argonne  National  Library  known  as 
LINPACK.  LINPACK  has  come  into  being  in  the  same  manner 
that  FUNPACK  did  cuid  is  well  documented  [36]  .  The  refer¬ 
enced  documentation  includes  a  listing  of  all  subroutines 
and  tables  of  timing  data  on  each  routine  for  each  of  the 
test  sites  [36 :B-1,D-11] . 

The  best  choice  of  these  programs  is  not  known. 
Since  the  second  program  never  fully  worked  and  produced 
correct  answers,  this  question  had  not  operationally  been 
decided.  The  plan  was  to  run  several  jobs  through  each  of 
the  programs  and  to  see  which  came  out  best  in  terms  of 
accuracy,  speed,  and  ease  of  operation.  The  result  of  the 
test  would  decide  which  routine  to  use  for  the  final  runs. 
For  the  problem  here,  time  is  the  biggest  problem,  not 


memory . 


The  purpose  of  the  input  program  was  to: 

1.  Read  in  the  geometrical  data  on  the  scatterer. 

2.  Calculate  the  elliptic  coordinate  system 
pareuneters . 

3.  Determine  the  endpoints  of  the  segments. 

4 .  Calculate  constemts  used  in  subsequent  pro¬ 
grams  . 

The  listing  is  given  in  Appendix  C.  The  calcula¬ 
tion  of  elliptic  coordinates  is  done  using  the  relation¬ 
ships  given  in  Appendix  A.  In  calculating  the  coordi¬ 
nates  for  the  line  source  location,  it  was  assumed  that 

the  sources  have  the  same  value  as  the  scatterer  does. 

o 

Therefore,  only  the  source  focal  length,  C_,  and  angle 

S 

coordinate,  v  ,  need  to  be  determined, 
s 

Wave  Number.  As  one  of  the  constants,  the  wave 
number  in  the  dielectric  was  determined.  It  has  a  value 
between,  and  k^,  the  wave  numbers  in  free  space  and  for 


Assuming  that  the  scatterer  will  locally  act  as  a 
flat  slab  waveguide,  the  electric  field  within  the  scat¬ 
terer  is  of  the  odd  type.  This  is  based  on  the  assumption 
that  the  electric  field  is  constant  through  the  dielectric 
shell  and,  therefore,  the  electric  field  is  symmetric 
about  the  middle  of  the  slab.  Since  there  is  no  variation 
in  the  electric  field  in  the  y  direction,  the  modes  excited 


vAiere 


(4.8) 


(4.9 


This  mecuis  that  an  approximate  value  for  k^j^T/2  is 


with  T  =  X  /20  and  e  =4 
o  r 

kj^T-  .5262 
k^X  =  10.5246 

(kyj^X)  2  =  ej.(k^X) 2  -  (k^X) 2  (4.1 

or 

kyj^X  *6.87  [37:12-13;15] 

The  determination  of  segment  and  points  was  a 
critical  phase  of  the  program.  No  segment  was  to  have  a 
length  greater  than  X/4  [29:5].  It  was  decided  that  for 
a  radius  of  curvature  of  2.5  or  greater,  that  this  segment 
length  was  adequate  since  a  straight  line  X/4  wavelength 
long  would  vary  little  from  the  arc  of  the  ellipse.  For 
the  segments  where  the  scatterer  has  a  smaller  radius  of 
curvature,  the  segments  were  to  be  p/10  wavelengths  long 
where  p  is  the  radius  of  curvature.  This  was  based  on  a 
straight  line  extrapolation  from  the  cutoff  point  of  .25 
for  2.5  (see  Figure  18). 

Based  on  this  the  cell  end  point  (one  was  known) 
was  determined  by  taking  the  arcsine  of  the  cell  length 
over  the  distance  from  the  origin  so  that  point  as  shown 


Figure  19.  Cell  Endpoint  Determination 

Matrix  Generation  and 
Algebra  Solver 

The  second  program  has  three  main  functions: 

1.  Determine  the  matrix  elements  of  the  reaction 
matrix,  Z^, 

2.  Determine  the  Vector  elements  for  Vu... 

3.  Solve  the  linear  algebra  problem, 

The  double  integral  generates  a  full  matrix  of  com¬ 
plex  elements  since  the  integral  over  the  nth  cell  (primed 
coordinates)  is  independent  of  the  integration  over  the 
mth  cell.  The  total  E  field  integral  generates  real  numbers 
only  and  contributes  only  to  the  matrix  diagonal  and  to  the 
major  codiagonals.  Due  to  the  E  field  integral,  a  blocked 
IF-THEM-ELSE  statement  is  used  to  test  which  contribution 


to  include.  The  program  calculates  the  real  and  the 
imaginary  parts  separately  because  the  IMSL  integrators 
are  real  functions. 


Algorithm.  The  sinusoidal  basis  function  had  to 

be  divided  up  into  two  subintervals;  v  .  to  and  to 

n— J.  n  n 

[24:5351  (see  also  Chapter  III).  In  the  double  inte¬ 
gral,  the  multiplication  of  the  sinusoid  for  the  mth  cell 
auid  the  sinusoid  for  the  nth  cell  results  in  four  double 
integrals  for  each  cell. 


^sin[k(v-v^j^)] 


sin[k(v^^-v)]  I 
+  sin[k(v^l-v)]  ( 


sintk(v'-v^_^)] 

*  sin(k(Vj^-v^^j^)  J 


+ 


sin[k(v^j^-v')]  j 
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V  V 
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sin  (k  ( v)  j  sin  [k  ( v ' )  ] 
sin(WA^j^-vJ  ]  sin[k  1  ^  ^ 


(4.11b) 


This  also  occurs  in  the  total  field  expansion  except  the 
integration  is  only  over  the  mth  cell  as  shovm  in  (4.11c) 


sin[k(v-Vj^j^)]  r  sin[k(v-v^_j^)]  sin[k(Vj^j^-v)]  ) 

''m-l 

^nH-1  sinlk(v^j^-v)]  <  sin[k(v-v^_j^)]  sin[k(Vj^j^-v)]  i 


(4.11c) 


m 


vdiere 

dA'  *  c^(cofih^  VI  -  cos^  V)  dV 
o 

2  2  2 
dA  =  c  (cosh  p  -  cos  v)  dv 
o 

Endpoints.  Figure  20  is  a  graphical  presentation 
of  the  double  integration  over  v  and  v' .  The  v  =  v'  line 
is  the  line  singularity  where  the  integrand  is  singular 
due  to  the  Bessel  function  of  the  second  kind.  Each  cell 
has  four  contributions  according  to  equation  (4.11b). 

This  can  be  seen  in  Figure  20.  As  an  example,  let  m  =  4 
and  n  =  2.  The  coverage  goes  back  to  m  =  3  and  up  to  m  =  5. 
The  integration  over  has  the  same  feature.  The  shaded 
block  represents  the  total  contribution  to  the  (4,2)  cell. 

Integration  on  the  cells  near  the  end  points 
presents  a  special  problem.  Since  the  scatterer  is  a  con¬ 
tinuous  object,  the  contribution  of  all  cells  must  be  the 
same.  The  affect  of  the  current  in  an  adjacent  region, 
but  beyond  the  symmetric  boundary  must  be  included. 
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Figure  21.  Current  Contribution  from 
Beyond  the  Line  of  Symmetry 

For  example/  the  (0/0)  cell  has  three  contributions  below 
the  eocis  in  Figure  20.  Originally/  the  program  handled 
each  case  on  the  end  points  as  a  separate  entity/  multiply¬ 
ing  the  (0/0)  and  (N/N)  cell  contribution  by  four  and  the 
ones  in  between  by  two  (point  G  on  Figure  20) .  However/ 
this  resulted  in  a  long/  bulky  program  and  it  was  difficult 
to  ascertain  if  the  contributions  from  cells  beyond  the 
symmetric  line  were  equal.  Thus  a  question  came  up  if 
multiplying  by  a  constant  would  produce  correct  results. 
Therefore/  a  "negative  cell"  was  created  by  folding  the 
cell  on  the  +/-90'’  line  over  the  other  side.  In  other 
words /  one  cell  was  created  that  was  equal  in  length  to 
the  0th  and  Nth  cell/  but  on  the  other  side  of  the  ±  90" 
lines  (see  Figure  21). 

To  assume  that  any  cells  on  the  other  side  of  the 


line  were  equal  in  length  to  the  adjacent  cell  was  an 


excellent  assuinption  since  the  region  of  concern  is  the 
essentially  flat  region  of  the  ellipse.  In  the  case  of  a 
circular  shell,  as  the  scatterer  becomes  more  circular, 
the  length  of  the  cells  becomes  more  nearly  equal.  This 
greatly  simplified  the  program  from  now  only  contribution 
from  the  total  field  expansion  had  to  be  checked. 

Additional  contributions  from  the  other  side  of 
the  ellipse  had  to  be  considered  due  to  the  assumed  sym¬ 
metry  of  the  problem  and  restriction  of  v  to  +/-90®  (see 
Figure  22).  Since  the  line  source  is  on  the  y  axis  (see 
Figure  6),  the  incident  energy  on  the  right  side  of  the 
shell  (+90  to  -90°)  is  equal  to  that  which  is  incident  on 
the  left  side  (180°  to  270°).  Therefore,  only  one-half  of 
the  ellipse  had  to  be  directly  evaluated,  thus  reducing 
the  number  of  cells  needed  by  a  half.  It  would  not  be 
correct  to  either  neglect  the  other  side  or  to  simply 
multiply  the  field  by  two. 

In  evaluating  the  Richmond  integral  equation 
(equation  (2.8)),  the  interaction  between  the  field 
generated  by  the  source  cell  and  the  observation  cell  is 
being  calculated  (see  Figure  8) .  This  must  also  be  con¬ 
sidered  with  the  image  cells  as  shown  in  Figure  22(a). 

This  was  done  in  the  program  by  integrating  over  a  cell 
180°  opposite  the  observation  cell. 

To  multiply  the  field  calculated  by  two  to  account 
for  the  image  contribution  would  give  erroneous  results 
for  the  distance  changes  as  the  position  on  the  ellipse 


changes  (see  Figure  22 (b)^.  The  use  of  symmetry  does  not 
reduce  the  number  of  calculations.  It  reduces  the  amount 
of  core  required  to  solve  the  problem/  a  critical  factor 
vdien  large  (greater  than  one  wavelength)  is  considered. 

Singularity.  The  progreun  had  to  evaluate  the 
principal  value  integral  around  the  singular  points.  A 
numerical  integration  routine  can  not  take  the  limiting 
value  of  a  function  and  special  routines  must  account  for 
these  regions.  In  this  problem,  ther^i  were  three  types 
of  singularities  encountered  ^ich  could  be  considered  as 
two  classes.  Table  5  gives  the  class  and  type  of  singular! 
ties  encountered.  Figure  23  is  a  graphic  presentation  of  . 
each  class. 


Table  5.  Program  Singularities 


Type  of 

Singularity 

Class 

Cross- 

Reference  to 
Figure  20 

Overlapping 

cells 

n  =  m 

line 

Region  B 

m  =  n  ±  1 

line 

Region  C 

Corner 

V  =  V 

m  n 

point 

Point  A 

There  were  three 

separate  methods 

tested  in  con- 

sideration  of  the  singular  cells  with  a  goal  that  accurate 
results  were  returned,  but  with  the  CPU  time  being  held  to 
a  minimtim.  The  first  routine  integrated  up  to  the  line 


singularity,  but  left  an  area  surrounding  the  singularity. 
Jeunes,  et  al.  [38:339-342]  suggested  that  the  function  be 


expanded  in  a  Taylor  Series  around  the  singular  point, 
integrating  term  by  term.  A  test  routine  was  written 
integrating  In  lx-x'|  over  -1  to  1  using  a  Sin?>son  Rule 
routine  [38:331]  and  the  IMSL  DCADRE  routine,  but  with  the 
singular  region  ignored.  Testing  only  10  subdivisions, 
the  error  was  -0.193  percent.  Using  1000  subdivisions, 
the  error  was  -2.5  x  lO”^  percent. 

The  point  singularity  was  handled  by  excluding 
the  point.  This  problem  occurred  when  the  range  of  each 
integration  overlapped,  say  going  from  1  to  2  and  2  to  3. 
Using  the  In  |x-x'|  function  again,  it  was  found  that  by 
integrating  the  inner  integral  from  1  to  1.999999999999 

Up 

and  then  the  outer  from  2  to  3,  the  results  were  limited 
only  by  the  limitation  on  the  accuracy  of  the  Simpson 
routine  and  the  ntimber  of  subdivisions  used. 

The  disadvantage  of  this  is  now  special  cases 
have  to  be  drawn  up  for  the  possible  singularities.  The 
comer  singularities  were  generally  found  by  trial  and 
error  (i.e.,  the  program  hit  a  divide  by  zero  error  and 
quit)  . 

The  first  method  discussed  involved  special  pro¬ 
grams  and  calls  to  avoid  the  problem.  In  either  case,  the 
singular  point  was  approached,  but  not  integrated  over. 
However,  it  was  decided  that  these  special  routines  were 


function  went  to  zero  in  the  program,  it  was  declared  to 

have  a  small  value.  The  minimum  value  that  the  FUNPACK 

—128  —39 

routine  could  handle  was  2  or  2.938  x  10 

-4  -25 

The  range  of  values  tried  was  from  10  to  10 
This  routine  was  very  easy  to  prograun;  however,  the  run 
times  doubled.  The  OCAORE  integrator  would  continuously 
subdivide  the  interval  from  the  last  actual  value  to  the 
declared  value.  Using  a  larger  number  did  not  resolve 
the  problem  since  the  integrator  would  bog  down  in  try¬ 
ing  to  resolve  the  sharp  cutoff  that  resulted. 

The  final  routine  involved  "capping"  the  singular 
point  with  a  parabolic  approximation.  This  cut  run  times 
by  one-half  since  now  the  singular  region  had  a  smooth 
peak,  as  opposed  to  a  sharp  peak  or’  flat  top. 

Equation  (4.12)  is  the  function  used.  Let 

H*(6)  =  a„  and  H*'(5)  =  a, 
o  1 


(4.12) 

The  (£)  used  was  10~^.  Using  a^  and  a^  with  (4.12)  results 
in 


(4.13) 


Using  this  function  sin^lified  the  problem  since  the 
result  could  be  used  by  either  class  of  singular  re^^ion. 
The  routine  significantly  reduced  run  times  also. 

Test  Runs.  To  test  the  program,  the  scatterer 
was  made  to  be  as  circular  as  possible.  The  resultamt 


6 


should  be  a  toeplitz  matrix.  In  such  a  matrix,  all 
of  the  elements  along  each  diagonal  are  equal.  The  toe¬ 
plitz  matrix  results  because  the  array  elements  are  cal¬ 
culated  on  a  basis  of  geometry  [14:341] .  All  of  these 
elements  are  the  result  of  integration  of  equally  sized 
self  cells.  Schneider's  results  showed  that  this  matrix 
does  result  [4:18].  The  test  program  included  a  repeat 
of  Schneider's  results  [4:27-30,47-49]  for  the  small 
scatterer.  This  would  validate  the  prograun  since  his 
results  are  correct  and  the  solution  is  unique  [40:532-534] 

Unresolved  Problems .  The  test  case  run  did  show 
that  the  toeplitz  matrix  does  result  for  the  circular 
scatter.  This  validates  the  fact  that  the  geometry  of 
the  problem  is  being  described  correctly  by  the  program. 

It  does  not  validate  that  the  values  obtained  are  correct. 
Herein  lies  the  problem.  The  resultant  reaction  matrix 
should  be  banded  in  that  the  elements  of  the  main  diagonal 
and  the  right  and  left  codiagonals  are  significantly  higher 
than  the  rest  of  the  array.  Additionally,  the  magnitude 
of  the  voltage  vector  should  increase  as  you  evaluate  from 
the  0th  cell  to  the  Nth  cell.  This  increase  should  be 
significant. 

Neither  of  the  described  conditions  was  met  by 
the  resultant  array  or  vector.  The  voltage  vector  gradu¬ 
ally  increases  as  you  get  closer  to  the  line  source,  but 
only  by  a  factor  of  3.  '"he  main  diagonals  are  only  an 
order  of  10  higher  thaui  the  other  elements  of  the  reaction 


array.  The  other  elements  then  varied  in  sign  and  magni> 
tude  across  a  row.  Their  magnitude  did  not  go  to  zero  as 
expected. 

Radiation  Calculation 

The  final  program  was  to  calculate  the  far  field 
electric  field  and  display  the  results.  The  current  in 
each  cell  would  be  obtained  by  going  back  to  the  expansion 
function  and  putting  the  calculated  back  into  the  equa¬ 
tion.  However,  the  far  field  is  very  insensitive  to  small 
changes  over  a  small  area.  Therefore,  the  electric  field 
for  the  cells  was  assvimed  to  be  V_,  (M)  across  the  entire 

MM 

mth  cell.  The  far  field  approximation  for  the  vector 
potential  is 


J  (p’)eSw’  eft'  (4.14) 

z 

[10:229] 

or  in  elliptic  coordinates,  if  p  is  the  distance  from  the 
origin  to  the  point  on  the  ellipse 


,2  -jkp  /ntl 


y  8  jifcp 


f  j  (y  ,v')e^^  [a  cos  (J>  +  b  sin  4)]^cos(v-v') 

i  A  II 


n 


d^ 
(4.15) 


a  s  semi-major  axis 


b  ■  semi-minor  axis 


»  tan  ( (a/b)  tanv) 

(fr*  =  tan~^  ((a/b)  tan  v*) 

2  2 

dA'  ®  (cosh  -  cos  V*)  dv ' 

2 

The  far  field  observation  point  is  2a  /X  [14:24].  The  far 
field  electric  field  is  then  determined  from 


or 


E,  =  juU_A 

Z  *^0  2 


[14:25] 


(4.16) 


2  4  3  V 

31c  c  oju  T-’e  °  r  nt  1 

E  =  - - /  J 

Z  _  /  I  Z  O 


V  8j7tk  P 


n 


.  ^jlc[a^cos^(|)  +  b^sin^4>]^  cos  (v-V) 


2  2 

(cosh  -  cos  V)  dv* 


(4.17) 


Since  program  2  never  worlced,  the  radiation  program 
was  not  developed.  The  contribution  from  each  cell  would 
be  summed  over  the  entire  360** .  The  plot  was  to  be  dis¬ 
played  on  the  off  line  calcomp  plotter  using  the  DISSPLA 
package . 


V.  Lessons  Learned 

The  exact  reason  that  the  program  failed  is  not 
known.  Sources  of  error  could  have  come  from  a  myriad  of 
sources.  Countless  checks  have  been  performed  on  the  pro¬ 
gram  and  no  error  can  be  found.  A  better  method  needs  to 
be  used  in  the  integration  if  the  method  of  moments  has  to 
be  used  to  solve  this  problem.  Blue's  article  showed  that 
the  number  of  cells  could  be  kept  small  by  using  different 
basis  functions  [26] .  A  fewer  number  of  cells  would  have 
substantially  improved  this  procedure  since  fewer  integra¬ 
tions  would  have  been  necessary. 

Despite  the  program  failure,  this  thesis  showed 
that  straight  application  of  the  moment  method,  even  with 
the  use  of  the  more  complicated  basis  functions,  is  not 
practical  for  larger  scatterers.  While  the  core  problem 
has  been  solved,  a  real  time  problem  has  surfaced.  If  it 
takes  1300  seconds  to  process  a  33  by  33  array,  a  120  x  120 
array  would  take  way  too  long.  The  expense  in  that  much 
conqputing  would  be  high. 

When  computing  with  wire  scatterers,  the  observa¬ 
tion  point  is  done  on  the  inside  of  the  metal,  while  the 
source  is  considered  on  the  outside  skin.  With  this  in 
mind,  the  observation  point  could  have  been  made  just 
above  u..  Since  the  electric  field,  and  therefore  the 


current/  is  constant  throughout,  the  electric  field  would 
not  have  changed.  The  runs  done  so  far  show  that  the 
singular  cells  take  10  seconds,  nonsingular  take  1.5  to 
0.5  seconds,  and  the  voltage  matrix  takes  less  than  0.1 
seconds.  Therefore,  the  problem  is  in  the  singularity. 

The  removal  of  the  singular  point  would  greatly  speed  opera 
tions.  While  time  would  not  permit  this  condition,  it  is 
worth  continued  investigation. 


VI.  Conclusions 


This  study  analyzed  the  scattering  of  cylindrical 
electromagnetic  waves  off  of  dielectric  scatterers.  The 
scatterer  was  an  elliptic  shell,  designed  to  model  the 
radome  of  the  F-16.  The  equation  developed  by  Richmond 
[10]  was  used  and  solved  by  the  method  of  moments.  To 
reduce  the  amount  of  storage  needed,  the  basis  and  test¬ 
ing  functions  were  the  sinusoidal  basis  functions.  The 
result  was  a  complicated  integral  that  took  a  great  deal 
of  time  to  compute  over  each  cell. 

The  progreun  never  worked  and  there  were  no  plots 
produced.  The  exact  source  of  error  is  unknown,  but  the 
most  probable  are  either  a  prograunming  or  an  error  in  the 
integration  over  the  shell  thickness.  The  study  did  show 
that  this  method  was  in^ractical  for  large  scatterers. 

The  eunount  of  integration  required  to  fill  the  reaction 
matrix  was  far  too  much  to  be  practical.  The  cost  of  such 
runs  made  justification  difficult. 

While  the  progreun  never  produced  valid  results, 
it  is  felt  that  the  conclusion  that  the  method  is  impracti 
cal  is  valid.  The  evaluation  of  the  individual  contribu¬ 
tions  would  still  require  16  integrations  to  be  done  per 
cell.  Each  integration  takes  considerable  amount  of  time. 
The  improbable  values  of  the  array  and  vector  elements  are 


more  likely  due  to  the  handling  of  the  singular  region 
rather  than  in  a  programming  error.  The  removal  of  the 
singularity  may  have  removed  the  main  reason  for  the  side 
lobe.  Thus  the  array  elements  were  only  gradually  varying. 

As  recommendations  for  future  study.  Blue's  method 
of  different  basis  fvinctions  bears  more  study.  The  reduc¬ 
tion  in  the  number  of  basis  functions  would  greatly  reduce 
the  number  of  integrations  needed. 

The  method  used  here  may  be  improved  by  breaking 
up  the  singular  cells  into  subregions.  If  the  singular 
area  was  square,  Richmond's  analytical  integration  [10:336] 
could  be  applied.  The  rest  of  the  cell  would  be  considered 
using  the  methods  described  in  this  report.  The  contribu¬ 
tions  would  then  be  added  together. 

Looking  at  other  methods,  if  an  efficient  program 
could  be  written  for  the  generation  of  Mathieu  functions, 
then  a  series  solution  would  be  possible.  Finally,  an 
asymptotic  evaluation  of  the  integrals  would  then  enable 
a  solution  to  the  integral  equation  to  be  computed  since 
the  integration  would  no  longer  be  necessary . 
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Appendix  A 

The  Elliptic  Coordinate  System 

The  elliptic-hyperbolic  coordinate  system  is  one 
of  the  eleven  orthogonal  coordinate  systems  which  is 
formed  from  first-  and  second-degree  surfaces.  This  appen 
dix  will  provide  a  con^endium  of  information  and  relation¬ 
ships  useful  in  analyzing  the  elliptic  shell.  The  best 
overall  source  for  coordinate  system  information  was 
Moon  and  Spencer's  Field  Theory  Handbook  [42] .  The 
Schaum's  Outline  Series  handbook.  Vector  Analysis  by 
Spiegel  [43]  and  Morse  and  Feshback's  Methods  of  Theo¬ 
retical  Physics  also  provide  some  valuable  insight  into 
the  system.  Moon  and  Spencer  and  Morse  and  Feshback  have 
considerable  discussion  on  the  separation  of  variables, 
especially  in  terms  of  the  Laplace  and  Helmholtz  equations 
Burnside  [44],  had  a  very  complete  list  of  relationships 
related  to  the  elliptical  geometry. 

Figure  A-1  shows  the  coordinate  system  relative 
to  the  xy  plane.  The  positive  z  axis  is  up,  out  of  the 
page.  The  defining  relationships  are 

X  a  c  cosh  ]i  cos  V 


y  a  c  sinh  y  sin  v 


(A.l) 


■c  > 


The  relationship  between  and  v  is 


(P  =  tan’’^[(b/a)  tanv] 


(A. 3) 


This  is  based  on  the  fact  that  a  »  c  cosh  V  and  b  » 
c  sinh  li.  The  range  of  v  is  from  0  to  360®. 
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The  unit  vectors  and  are 
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2  2  2  2 
cosh  y  sin  v  +  sinh  y  cos  v 


^  is  the  vector  tangent  to  the  ellipse  and  y^  is  the 
vector  tangent  to  the  hyperbola  or  the  outward  normal  from 
the  ellipse  [44:310-311]. 

The  following  is  then  a  list  of  relationships  used 


or  noted  during  the  course  of  this  research: 


Differential  area 


P 


2  2  2 
dA  «  c  (cosh  y  -  cos  V  )  dydv 

2.2  2 
dA  =  c  (sinh  y  +  sin  v  )  dydv 

These  equations  are  equivalent. 

Elliptic  cylinders: 

2  2 
(x/c  cosh  y )  +  (y/c  sinh  y )  *  1 

arc  length: 


[42:18]  (A. 5a) 
[43:139]  (A. 5b) 

[42:17]  (A. 6) 
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Circumference : 

L  »  4aE(a/c) 

=  7r[3(a+b)  -  •y(a+3b)  (3a+b)] 
s  2TrVl/2(a^+b2) 
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[46:18]  (A. 8) 
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E(a/c)  is  the  complete  elliptic  integral  of.  the  second 
kind. 

Area 


A  =  irab  [46:18]  (A. 9) 

Radius  of  curvature: 

3/2 

p  -  (a^  sin^  V  +  b^  cos^v)  /ab  [46:21]  (A. 10) 
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See  [42]  for  equations  related  to  the  separation  of 
Laplace's  equation  and  the  Helmhotz  Equation  [42:18-20]. 
See  [45]  and  [46]  for  the  imsre  novel  relationships  for 
ellipse.  Lockwood  also  has  a  great  deal  of  infomation  on 
the  geometrical  properties  of  curves  in  general. 


Appendix  B 
FUNPACK  Release  2 

The  following  is  a  copy  of  a  listing  executed  from 
the  EDIT  LIB  Users  Library.  Since  the  special  functions 
were  all  generated  using  the  FUNPACK  library,  and  since 
there  is  not  a  commercially  published  manual  available, 
this  file  is  included  as  part  of  this  thesis.  More  detailed 
information  on  FUNPACK  is  availadile  in  AFFDL-TM-77-89-FBR 
[33].  For  detailed  data  on  LINPACK  and  the  IMSL  see  [36] 
and  [35] ,  respectively. 


LFUN 


THIS  LISTING  IS  OUTPUT  FROM  PROCEDURE  LFUN(LFN)  EXECUTED  FROM  THE 
EDITLIB  USER  LIBRARY.  IT  DOCUMENTS  COOES  SELECTED  FROM  ARCONNE  CODE 
CENTER  NO.  610,  FUNPACK  RELEASE  2,  TO  EVALUATE  CERTAIN  SPECIAL 
FUNCTIONS.  THE  CODES  DESCRIBED  ARE  IN  AN  EDITLIB  USER  LIBRARY  AS 
CENTRAL  PROCESSOR  PROGRAMS  COMPILED  UNDER  FORTRAN  EXTENDED,  VERSION 
4.5+414,  USING  THE  ROUNDED  ARITHMETIC  OPTION.  THE  SOURCE  CODES  FROM 
WHICH  THEY  ARE  DERIVED  ARE  IN  AN  UPDATE  OLDPL.  BOTH  THE  USER  LIBRARY 
AND  OLDPL  RESIDE  ON  MAGNETIC  TAPE.  FOR  FURTHER  INFORMATION  ON  THESE 
CODES  OR  TO  ACCESS  THEM  FROM  TAPE,  CALL 

DONALD  S.  CLEMM  /  AFWAL/FIBR  /  513-255-5350  (Av  735-5350/. 

********************************************************************** 

ACC  ABSTRACT  610 

1.  NAME  OR  DESIGNATION  OF  PROGRAM  -  FUNPACK  RELEASE  2 

2.  COMPUTER  FOR  WHICH  PROGRAM  IS  DESIGNED  AND  OTHERS  UPON  WHICH 
IT  IS  OPERABLE  -  IBM360,370,  CDC6000-7000,  UNIVAC1108, 1110 

3.  DESCRIPTION  OF  PROBLEM  OR  FUNCTION  -  FUNPACK  IS  A  COLLECTION  OF 
FORTRAN  SUBROUTINES  TO  EVALUATE  CERTAIN  SPECIAL  FUNCTIONS.  THE 
INDIVIDUAL  SUBROUTINES  ARE  - 


IDENTIFICATION 

DESCRIPTION 

NATSIO 

F2I0 

BESSEL  FUNCTION  I-SUB-0 

NATSIl 

F2I1 

BESSEL  FUNCTION  I-SUB-1 

NATSJO 

F2J0 

BESSEL  FUNCTION  J-SUB-0 

NATSJl 

F2J1 

BESSEL  FUNCTION  J-SUB-1 

NATSKO 

F2K0 

BESSEL  FUNCTION  K-SUB-O 

NATSKl 

F2K1 

BESSEL  FUNCTION  K-SUB-l 

NATSBESY 

F2BY 

BESSEL  FUNCTION  Y-SUB-NU 

DAW 

FIDW 

DAWSON'S  INTEGRAL 

ELIPK 

FIEK 

COMPLETE  ELLIPTIC  INTEGRAL 

OF 

THE 

FIRST  KIND 

ELIPE 

FIEE 

COMPLETE  ELLIPTIC  INTEGRAL 

OF 

THE 

SECOND  KINO 

El 

FIEI 

EXPONENTIAL  INTEGRALS 

NATSPSI 

F2PS 

PSI  (LOGARITHMIC  DERIVATIVE  OF  GAMMA  FUNCTION) 

MONERR 

FIMO 

ERROR  MONITORING  PACKAGE 

4.  METHOD  OF  SOLUTION  -  FUNPACK  USES  EVALUATION  OF  MINIMAX  APPROXI¬ 
MATIONS. 

5.  RESTRICTIONS  ON  THE  COMPLEXITY  OF  THE  PROBLEM  - 

6.  TYPICAL  RUNNING  TLME  - 

7.  UNUSUAL  FEATURES  OF  THE  PROGRAM  -  THESE  ROUTINES  HAVE  BEEN 
CERTIFIED  UNDER  THE  NATS  PROJECT  FOR  THE  MACHINES  AND  OPERATING 
SYSTEMS  INDICATED  IN  ITEM  13  AND  FOR  THE  COMPILERS  INDICATED  IN 
ITEM  12.  EXTENSIVE  TESTING  ON  THESE  MACHINES  HAS  SHOWN  NO  EVI¬ 
DENCE  OF  PERFORMANCE  DIFFICULTIES.  EXCEPTIONS,  IF  ANY,  FOLLOW  - 

CDC  VERSIONS  OF  THESE  SUBROUTINES  ARE  TUNED  TO  PERFORM  BEST 
USING  THE  ROUNDED  ARITilMETIC  OPTION  ON  CDC  COMPILERS. 

THE  ACCURACY  OF  THE  SUBROUTINES  FOR  THE  ELEMENTARY  FUNCTIONS 
(EXP,  ALOG,  ETC.)  CAN  AFFECT  THE  ACCURACY  OF  FUNPACK  SUBROUTINES. 

THE  IBM  VERSION  OF  THIS  PACKAGE  ASSUMES  THE  IBM-SUPPLIED 
TRACEBACK  SUBROUTINE  ERRTRA  IS  AVAILABLE. 

THE  NATS  PROJECT  FULLY  SUPPORTS  CERTIFIED  ROUTINES  IN  THE 
SENSE  THAT  REPORTS  OF  POOR  OR  INCORRECT  PERFORMANCE  ON  AT  LEAST 
THE  MACHINES  AND  OPERATING  SYSTEMS  LISTED  WILL  BE  EXAMINED  AND 
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8. 


9. 


10. 


11. 

12. 


NECESSARY  CORRECTIONS  MADE.  THIS  ASSURANCE  OF  SUPPORT  APPLIES 
ONLY  WHEN  THE  SOFTWARE  IS  OBTAINED  DIRECTLY  FROM  THE  AKCONNE 
C01»  CENTER  AND  HAS  NOT  BEEN  MODIFIED. 

REUTED  AND  AUXILIARY  PROGRAMS  -  FUNPACK  RELEASE  2  REPLACES 
FUNPACK,  AN  EARLIER  PACKAGE  SUBMITTED  IN  JULY  1973  AND  DISTRI¬ 
BUTED  BY  THE  CODE  CENTER  AS  ACC  NO.  610  PRIOR  TO  THIS  RELEASE. 
SUBROUTINES  IDENTIFIED  AS  FIXX  IN  ITEM  3  ABOVE  ARE  UNMODIFIED 
FROM  THE  PREVIOUS  RELEASE.  WITH  THE  POSSIBLE  EXCEPTION  OF  THE 
TEST  MATERIAL  FOR  THE  IBM  MACHINES. 

STATUS  -  ABSTRACT  FIRST  DISTRIBUTED  JULY  1973. 

IBM360.370  VERSION  OF  FUNPACK  SUBMITTED  AUGUST  1973, 

■  REPLACED  BY  FUNPACK  RELEASE  2  SEPTEMBER  1976,  SAMPLE 
PROBLEMS  EXECUTED  BY  ACC  SEPTEMBER  1976  ON  AN 
IBM370/195. 

COC6000-7000  VERSION  OF  FUNPACK  SUBMITTED  AUGUST  1973, 
REPLACED  BY  FUNPACK  RELEASE  2  SEPTEMBER  1976,  SAMPLE 
PROBLEMS  EXECUTED  BY  ACC  OCTOBER  1976. 

UNIVAC1108  VERSION  OF  FUNPACK  SUBMITTED  AUGUST  1973, 
REPLACED  BY  UNIVAC1108, 1110  VERSION  OF  FUNPACK 
RELEASE  2  SEPTEMBER  1976. 

REFERENCES  -  J.  M.  BLAIR  AND  C.  A.  EDWARDS,  STABLE  RATIONAL 
MINIMAX  APPROXIMATIONS  TO  THE  MODIFIED  BESSEL  FUNCTIONS  I0(X)  AND 
I1(X),  AECL-4928,  1974. 

J.  M.  BLAIR  AND  A.  E.  RUSSON,  RATIONAL  FUNCTION 
MINIMAX  APPROXIMATIONS  FOR  THE  BESSEL  FUNCTIONS  KO(X)  AND  K1(X), 
AECL-3461,  1969. 

W.  J.  CODY,  CHEBYSHEV  APPROXIMATIONS  FOR  THE  COM¬ 
PLETE  ELLIPTIC  INTEGRALS  K  AND  E,  MATH.  COMP.  19,  105-112  (1965). 

W.  J.  CODY,  R.  M.  MOTLEY,  AND  L.  W.  FULLERTON,  THE 
COMPUTATION  OF  REAL  FRACTIONAL  ORDER  BESSEL  FUNCTIONS  OF  THE 


SECOND  KIND,  APPLIED  MATHEMATICS  DIVISION  TECHNICAL  MEMORANDUM 
NO.  291,  ANL,  1976. 

W.  J.  CODY,  K.  A.  PACIOREK,  AND  H.  C.  THACHER,  JR., 
CHEBYSHEV  APPROXIMATIONS  FOR  DiMSON'S  INTEGRAL,  MATH.  COMP.  24, 
171-178  (1970). 

W.  J.  CODY,  A.  J.  STRECOK,  AND  H.  C.  THACHER,  JR., 
CHEBYSHEV  APPROXIMATIONS  FOR  THE  PSI  FUNCTION,  MATH.  COMP.  27, 
123-127  (1973). 


W.  J.  CODY  AND  HENRY  C.  THACHER,  JR.,  RATIONAL 
CHEBYSHEV  APPROXIMATIONS  FOR  THE  EXPONENTIAL  INTEGRAL  E1(X), 
MATH.  COMP.  22,  641-649  (1968). 

W.  J.  CODY  AND  HENRY  C.  THACHER,  JR.,  RATIONAL 
CHEBYSHEV  APPROXIMATIONS  FOR  THE  EXPONENTIAL  INTEGRAL  EI(X), 


MATH.  COMP.  23,  289-303  (1969). 
MACHINE  REQUIREMENTS  - 
PROGRAMMING  LANGUAGES  USED  - 
FORTRAN  IV(C) 

FORTRAN  IV(G  21) 

FORTRAN  IV(Gl) 

FORTRAN  IV(H) 

FORTRAN  IV(H  20.1) 

FORTRAN  IV(H  21.6} 

FORTRAN  IV(H  21.7) 

FORTRAN  IV(U  21.8) 


IBM360/75, 195, IBM370/165, 195 
IBM360/6  5,75,91, IBM3  70/158 
IBM360/75,195 

IBM360/65,67, 75, 195, IBM370/165 

IBM360/6S, 75,195,AM0AHL470V/6 

IflM360/75,91 

IBM360/75 

IBM360/6S 
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FORTRAN  1V(H  EXTENDED) 
FORTRAN  IV(H  EXTENDED  2.1) 
FORTRAN  IV(WATFIV) 

FORTRAN  IV(UATF1V  1.4) 

FTN 

FTN(4.642) 

FTN(3.0) 

RUN 

RUN(2.3) 

FUN 
FTN  V 

FTN  V(SIQA-O) 

FTN  V(MACC  1.175) 


1BM370/195 

1BM360/75, 91, IBM370/168, 168-11 

IBM360/67 

IBM370/158 

CDC6400.6S00 

CDC6600,CY173,175 

CDC6400 

CDC6400. 6500, 6600-6400, 7600 
CDC6400 

CDC6400-6500, 6600, 6700 
UNIVAC1108 
UNIVAC1108 
UNlVAClltO 


OPERATING  SYSTEM  OR  MONITOR  UNDER  WHICH  PROGRAM  IS  EXECUTED  - 


OS/360 

OS/360(19.6) 

OS/360(20.1) 

OS/360(20.7) 

OS/360(21.0) 

0S/360(21.7) 

OS/360(21.8) 

OS/MVT(21.7) 

OS/MVT(21.8) 

0S/VS2(1.6) 

MTS 

STANFORD  UNIVERSITY 
PURDUE  UNIVERSITY 
BERKELEY  LABORATORY 
LIVERMORE  LABORATORY 
NCAR 

SCOPE(3.3) 

SCOPE(3.4) 

UT2D 
EXEC  8 
EXEC  31.244E 
EXEC  MACC  31.66 


IBM360/67 

IBM360/65 

IBM370/165 

IBM360/75,195 

1BM360/75,91 

IBM360/75, 370/195 

IBM360/75 

IBM370/165-I1 

IBM360/6S, 91 , IBM370/I 58 

IBM370/168 

IBM360/67,AMDAHL470V/6 

IBM360/67 

CDC6500, 6400-6500 

CDC6600,7600 

CDC6600,7600 

CDC6600, 7600 

CDC6400 

CDC6600 

CDC6600-6400 

UNIVAC1108 

UNIVAC1108 

UNIVAClllO 


OTHER  PROGRAMMING  OR  OPERATING  INFORMATION  OR  RESTRICTIONS  - 
LOCATIONS  AND  MACHINES  USED  FOR  FUNPACK  TESTING  WERE  - 


MACHINE  TEST  SITE 

IBM360/65,IBM370/158  AMES  LABORATORY,  IOWA  STATE 

UNIVERSITY 

IBM360/75, 195, 370/195  ARCONNE  NATIONAL  LABORATORY 

1BM360/75,91  OAK  RIDGE  tMTIONAL  LABORATORY 

IBN360/67,91,370/168  STANFORD  UNIVERSITY 

IBH360/75  STOCKHOLM  DATA  CENTER 

1BM360/65  THE  UNIVERSITY  OF  CHICAGO 

IBH360/75  UNIVERSITY  OF  ILLINOIS  AT 

UKBANA-CHAMPAICN 

IBM360/67,AMDAUL470V/6  THE  UNIVERSITY  OF  MICHIGAN 

IBN360/67  THE  UNIVERSITY  OF  NEW  MEXICO 

IBM370/165, 165-11  UNIVERSITY  OF  TORONTO 

CDC6600  KIRTLAND  AIR  FORCE  RASE/AFWL 

CY173,175  ICASE/NASA  LANGLEY  RESEARCH  CENTER 

CDC6600,7600  LAWRENCE  BERKELEY  LABORATORY 


1BM360/75,91 

IBN360/67,91,370/168 

IBH360/75 

1BM360/65 

IBH360/75 
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CDC6600,7600  '  LAWRENCE  LIVERMORE  LABORATORY 

CDC6600,7600  NATIONAL  CENTER  FOR  ATMOSPHERIC 

RESEARCH 

CDC6400  NORTHWESTERN  UNIVERSITY 

CDC6A00-6500  PURDUE  UNIVERSITY 

CDC6600-6400  THE  UNIVERSITY  OF  TEXAS  AT  AUSTIN 

UNIVAC1108  ILLINOIS  INSTITUTE  OF  TECHNOLOGY 

UNIVAC1108  JET  PROPULSION  LABORATORY 

UNIVACl 108, 1110  UNIVERSITY  OF  WISCONSIN 

15.  NAME  AND  ESTABLISHMENT  OF  AUTHOR  - 

W.  J.  CDDY 

CONTACT  BURTON  S.  CARBOW 

APPLIED  MATHEMATICS  DIVISION 
ARGONNE  NATIONAL  LABORATORY 
9700  SOUTH  CASS  AVENUE 
ARGONNE,  ILLINOIS  60439 

16.  MATERIAL  AVAILABLE  >  MAGNETIC  TAPE  TRANSMITTAL 

SOURCE  DECKS  (370-3902  CiUlDS,  7600-3983  CARDS,  1110-411S 
CAROS) 

DEMONSTRATION  PROGRAM  SOURCE  DECKS  (370-2083  CAROS,  7600-2151 
CAROS,  1110-2211  CARDS) 

DEMONSTRATION  OUTPUT  (370-61  PAGES,  7600-58  PAGES,  1110-44 
PACES) 

MACHINE -READABLE  DOCUMENTATION  (370-4532  CARDS,  7600-4518 
CARDS,  1110-4518  CAROS) 

17.  CATEGORY  -  P 

KEYWORDS  -  SPECIAL  FUNCTIONS,  EXPONENTIAL  ItrTECRALS,  COMPLETE  ' 
ELLIPTIC  lOTEGRALS,  DAWSON'S  INTEGRAL,  BESSEL 
FUNCTIONS,  NEUMANN  FUNCTIONS,  PSI  FUNCTION 


Input  and  Segmentation  Program 


The  following  is  a  copy  of  the  source  code  for 
program  1.  Detailed  description  of  the  progreua  is  pro 
vided  in  Chapter  IV  of  this  report  and  in  the  comments 
located  within  the  code. 


APPENDIX  C.  COOHDIKATE  system  and  SEGMEirTATIION  PROGRAM  PACE  C-  1 


PROCRA? :  READI M  ( I MPUT ,  OUTPDT ,  DATA ,  TAPE  5-IMPUT ,  TAPE  6-OirrPUT , 

1  TAPE 7-DATA) 

Tins  SECTION  WILL  READ  IN  THE  VALL’ES  FOR  THE  SIZE  OF  THE 
ELLIPTIC  SHELL,  CatPUTE  THE  ELLIPTIC  COORDINATES,  OUTPUT 
THOSE  VALUES,  PLOT  THE  CYLINDER,  AND  S13:TC!1  IK  (WITH  THE 
SECEMENTS  KIPIBERED)  THE  LOCATION  OF  THE  SEGMENTS. 

DIMENSION  RAOIKO;  100), Ra^CO:  100),  PHI(0;100) 

REAL  K,KN0T,LA!1DA 

REAL  MUKOT,  MUIN,  IRJOL’T,  LENSIC(0: 100) ,  KUSEG(0:100) 

PI  -  3.1415926535898 

RL\D( 5 , 100 ,END-9999)A ,B ,T, XS , YS , FREQ 

FORMAT(6F10.6) 

IF  (A.EQ.E)  THEN 

B  -  A  *  0.999999999999 
E:n)  IF 

C  -  SQRT(A**2  -  B**2) 

E  -  C  /  A 

CIRCUM  -  4.0  *  A  *  ELIEl(E) 

DISTSC  -  S0RT(XS**2  +  YS**2) 

ICNOT  «  2.0  *  PI  *  FREQ  *  1.0E09  /  2.997925E03 
MUNOT  -  ATAJni(B/A) 

TAU  -  T/A 

IF  ((XS.EQ.O.).ANO.(YS.GT.O.))  TIIEH  . 

THETAS  -  90.0 
VS  -  90.0 

ELSE  IF  ((XS.EQ.O.).AKD.(YS.LT.O.))  THEN 
THETAS  -  270.0 
VS  -  270.0 

ELSE  IF  ((XS.EQ.O.).A.ND.(YS.EQ.O.))  THEN 
THETAS  -  0.0 
VS  -  0.0 
ELSE 

THETAS  -  ATA:;(YS/XS)  *  (180.0/PI) 

VS  •  atai:(a*tai:d(thetas)/b)  *  (iso.o/pi) 

END  IF 

IF  ((XS.Eq.O.).AJ;D.(YS.EQ.O.O))  THEN 
CS  -  0.0 

ELSE  IF  (XS.EO.O.O)  THEN 

CS  -  YS/(SIWKI!UN0T)*SIND(VS)) 

ELSE 

CS  -  XS/(COSH(MUNOT)*COSD(VS)) 

END  IF 

WAVELT  -  2.997925E08/(FPJ:Q*1.0E09) 

AW  -  A/WAVELT 
BW  -  B/WAVTLT 
CV7  -  C/WAVELT 
TW  -  T/ WAVELT 

K  -  SQRT(16.0*(PI**2)-(80.0*SQRT(3.0*(?I**2)+100.0)-800.0))/WAVELT 


APPENDIX  C.  COOl^niNATS  SYSTEI-!  AKD  SEGr!E:rrATI10N  PIIOGRAM  PACE  C- 


LAMDA  -  0.75  *  (V.*C)**2 

FC  -  K  *  C  *  Sl^^!(2.0*:lU^’0T)  *  (TAU  *  COSil(TAU)  -  Sm(TAU))  /S.O 
FCONSl  -  -(KKOT**2)/(8.0  *  PI  *  FREO  *  P.854185E-3) 
rCONS2  -  (KNOT**3)  *  (Sira(2.0  *  !!UKOT)**2) 

1  *  (TAU  *  COSII(TAU)  -  SIWKTAU)) 

2  /  <64.0  *  PI  *  FREQ  *  8.854185E-03) 

MUIN  -  MU^:OT  -  TAU/2.0 

HUOUT  -  :1U!:0T+  TAU/2.0 
MRITE(6,1030) 

1030  F0miAT(’l*,132('*')) 

URITE(6,1040)  A,  T,  AW,  T\^,  B,  EW,  XS,  YS,  C,  E,  CIRCUM 
1040  F0R:fAT('0',15X, 'INPUT  DATA: '/ESX, 'SEIII-MAJOR  AXIS:  ',n0.6, 

1*  ^^CTF.RS. '  ,10X, 'SHELL  THICra.’ESS:  ',F10.6,'  METERS.'/ 

A  30X,'(',F10.6,'  WAVELENGTHS) ',15X,'(’,F10. 6,'  WAVELENGTHS)'/ 
225X, 'SEMI-MINOR  AXIS:  ',F10.6,'  '1ETERS. ' , lOX, ' ( ' ,F10.6, '  WAVE-' 

B  , 'LENGTHS) '/25X, 'SOURCE  X', 

4'  COORDINATE:  ',F10.6,'  METERS. ' ,10X, 'SOURCE  Y  COORDINATE:  ', 
5F10.6,'  IIETEP.S.'//15X,'CALCULAT' 

6, 'ED  DATA: '/25X, 'FOCAL  LENGTH: ' ,5X,F10. 6, '  PIETERS. ', lOX, 

7 'ECCENTRICITY  :  '  ,F10.6,  *  . ' /25X, 'CIRCU^lFEPJcNCE: 

6F10.6,'  METERS.') 

WRITE(6, 1050)  )(UOUT,TAU,!tUIN,MUO’JT 
1050  F0RMAT('-',15X,•C0^’STA^T  ELLIPTIC  COORDINATES: '//,25X, 

I'ffUN'OT;  ',F10.6,'.',10X,’TAU(-T/A):  '  ,F1C. 5/2 5X, 'INNER  ', 

2 'RADIUS:  ' ,Fi0.6,10X, 'OUTER  RADIUS:  ',F10.6,'.') 

WRITE(6,1055)  ^^JOUT,  VS,  DISTSC,  THETAS,  CS 
1055  F0RMAT('0',15X, 'SOURCE  LOCATION: '/I 5X, ' (SECONT)  SOURCE  LOCATED  ' 
1,'180  DECREES  FROM  THE  ONE  INPUTTED) ' //25X, 'ELLIPTIC  COORDINATE:' 
2, 2X,F10. 6, lOX, 'ELLIPTIC  ANCL’LAR  COORDINATE:  ',rl0.6,'  DEGREES.'/ 
325X, 'POLAR  DISTANCE  FROT!  CENTER:  ',F10.6,'  METERS. ', lOX, 'POLAR  ' 
4'ANGLE:  ',F10.6,'  DEGREES. ' /50X, 'FOCAL  LENGTH  FOR  SOURCE  ELLIPSE 

5':  ',F10.6,'  METERS.'///) 

WRITE<6,1065)  FREO,WAVELT,K,Ki:OT 

1065  FOPJtATC  ',15X,'ELECTRaiAGKETIC  PAPJUIETERS'/25X, 'FREQUENCY:  ', 
1F10.6,'  GIGAHERTZ. ',1(K,' WAVELENGTH:  ',F10.6,'  METERS. ' /25X, 

2 'WAVE  NUMBER(OIELECTRIC):  ' ,F10.6, ' . ' ,10X, 'WAVE  OTMEERCFREE  ', 
3'SPACS):  ',F10.6, '.'////) 


CONCERT  ALL  METER  tlEASURES  TO  WAITLENGTH  MEASURES. 

CIRCUM  -  (CIRCUM/ 2. 0)/:;avelt 
J-0 

PHI(O)— 90.0 
LENSIG(0)-0.0 
NUSEG(O)— 90.0 
CURLEN-0.0 
THETA  -  -90.0 
50  COrJTIN’UE 

RO?J(  J  )-(  (  AW*  AW*S  IND(THETA  )  *S  INn(T:ir.TA  )  )+(  BW*BW*COSD  (THETA  )  * 


ArPEtJDIX  C.  COORDIKATC  SYSTEM  AJID  SEO'ECTATIIOR  PROGRAM  PAGE  C- 


1  C0SD(T11ETA)))**1.5/(AW*P.W) 

rw\DII(J)-SORT(AW*A»J*COSD(THETA)*COSO(THETA)+ 

1  BU*BU*Si:JD(TllETA)*SIMa{TllETA) ) 

J  ■  J  +  1 

ir(ROW(J-l).GE.2.5)  THEN 
LEMSIG(J)  -  .25 
ELSE 

LENSIC(J)  -  ROW(J-1)/10.0 
EfJD  IF 

C1!I-ATAN’(LEKSIG(J)7RADII(J-1  ))  *180. 0/PI 
CURLEN  •  CURLEK  +  LEKSIC(J) 

THETA-THETA+CIII 

PIII(J)-TIIETA 

IF  (PHI(J).EQ.-90.0)  THEN 
NUSEC(J)  -  PKI(J) 

ELSE  IF  (PHI(J).EQ.90.0)  THEN 
IRJSEGCJ)  -  PHI(J) 

ELSE 

KUSEC(J)  -  ATAN((A/B)*TA?:d(PHI(J)))  *  180.0/PI 
ENT)  IF 

IF  (THETA.lt. 90.)  THEN 
GO  TO  50 
ELSE 

LE:JSIG(J  )-CIRCU[l-CURLEN 
PllI(J)-90.0 
:juseg(j)-9o.o 
END  IF 

WRITE (6, 1060) 

DO  60  KK-0,J 

WRITE(6 , 1070)KK,LE>r,IG(KK) ,PKI(KK) ,NUSEG(RK) 

60  coirriNUE 

WRITE(6,1080) 

WRITE (6, 1090) 

1050  F0RMAT('0',5X, 'SEGMENT  rJUMBEPx '  ,6X, 'SEGMErT  LENGTH' ,8X, 'END  POI 

1,'KT  ANGLE ',1 IX, 'END  POINT  ANCLE ’/26X, ' (WAVELENGTHS) ' , 

26X, ' (DEGREES-CYLIMDRICAL) ' ,6X, ' (DEGRCES-ELLIPTICAL) ' //) 

1070  FORMATC  ' ,llX,I3,13X, F10.6, 12X, F10.6, 12X, F10.6) 

1080  FORMAT('l') 

1090  F0RMAT('+',25X,'END  OF  JOB.') 

WRITE ( 7  )  PI ,  riUKOT , TAU ,  VS ,  CS,  KNOT ,  K ,  WA VELT ,  LA'IDA ,  FC ,  FCOrS  1 ,  FCO.NS 2 , 
1  C,J,NUSEG,PHI,A,B,T,FREQ 
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EliCTRnriAC.MF.TIC  PARAHETERS 

FP.KQUF.rJCY;  2.000000  CICAIIERTZ.  WAVELENCTH:  .IA9«%  METERS. 

\JAVI.  NliriBER(l)IEI,ECTRTr);  AS. 807338.  WAVE  mMBER(FRE£  SPACE):  Al. 916895 


OT  IIWIRER 

SEGMENT  LENGTH 
(WU’ELENCTIIS) 

END  POINT  ANCLE 
( DEGREES-CYLINDRICAL) 

EN2)  POINT  ANCLE 
( DSGREES-ELLIPT IC AL ) 

0 

0.000000 

-90.000000 

-90.000000 

1 

.025000 

-84.289407 

-34.289407 

2 

.025000 

-78.578314 

-7S.578C14 

3 

.025000 

-72.868221 

-72.868221 

4 

.025000 

-67.157627 

-67.157627 

5 

.025000 

-61.447034 

-61.447034 

6 

.025000 

-55.736441 

-55.736441 

7 

.025000 

-50.025848 

-50.025848 

S 

.025000 

-44.315255 

-44.315255 

9 

.025000 

-38.604662 

-38.604662 

10 

.025000 

-32.894069 

-32.894069 

11 

.025000 

-27.183475 

-27.183475 

12 

.025000 

-21.472882 

-21.472882 

13 

.025000 

-15.762289 

-15.762289 

14 

.025000 

-10.051696 

-10.051696 

15 

.025000 

-4.341103 

-4.341103 

16 

.025000 

1.369490 

1.369490 

17 

.025000 

7.080083 

7.080083 

18 

.025000 

12.790676 

12.790676 

19 

.025000 

18.501270 

18.501270 

20 

.025000 

24.211863 

24.211863 

21 

.025000 

29.922456 

29.922456 

22 

.025000 

35.633049 

35.633049 

23 

.025000 

41.343542 

41.343642 

24 

.025000 

47.054235 

47.054235 

25 

.025000 

52.764828 

52.764828 

26 

.025000 

58.475422 

58.475422 

27 

.025000 

64.186015 

64.186015 

28 

.025000 

69.896608 

69.896608 

29 

.025000 

75.607201 

75.607201 

30 

.025000 

81.317794 

81.317794 

31 

.025000 

87.028387 

87.028337 

32 

-.014602 

90.000000 

90.000000 
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MATRIX  CALCULATOR  AMO  LI^X^R  ALGE»>RA  EQUATIOM  SOLVER  PACE  O' 


PROGRAM  ARRAY(DATA .OirrPUT, RESULT, TAPE 7-nATA,TAPC6-Ol!TPUT, 

1  TAPE8-RESULT) 

C 

C  THIS  SLTJROUTlfrES  CALCULATES  THE  FIELD  ARRAY,  ZMN,  AMD  THE 

C  SOURCE  VECTOR,  VMI,'.  THIS  IS  DOME  BY  FIRST  CALCULATING  THE 

C  ZMN  TERMS  WITH  M  BEING  HELD  CONSTANT.  EEFOPvE  LOOPING  TILE  M 

C  DO  LOOP,  THE  COORESPOMDIMG  VOLTAGE  VECTOR  VALUE  IS  CALCULATED 

C  BY  CALLING  THE  SUBROUTir:S  KNOWN  AS  “VOLTS". 

C 

CmMOK/CSLLS  /MU5 ,  NU6 ,  !IU1 ,  MU2 ,  NU3 ,  NU4 ,  M ,  N,  NOSEG ,  NUS  EC 
COMMON/ELLIPS/LPJ ,  C ,  TAU ,  K 

COMMOLVCONST/IFLAC, FC,FCONSl , FCOMS2, PI ,LAMDA 
COMMOM/SOURCE/CS,JO;OT,VS  ' 

CO:  1MON/  S  INGLE  /AZERO ,  AO^rE 

REAL  K,KNOT,LAMDA,MU,KUSEGD,t:USEC(0: 32) ,NU1  ,KU2,NU3,::U4,NU5,r:US 
CaiPLEX  ;  32 , 0: 32 ) , \R1M(0 : 32 ) 

DIMENSION  IPVT(32),Z(32) 

C 

C.\LL  UERSET(1,IX) 

READ(7 )  PI  ,*rU, TAU, VS ,  CS, IDfOT, K.UAVELT, LAMDA,  FC, FCONSl ,  FC0NS2 , 

1  C, NOS EG, NUS EG 

.PRINT  ’,PI,'  ’.TAU,’  ’,VS,’  ’,CS,’  ’,K,\’OT 

PRIOT  *,’  ’,K, '  ’,UAVELT,’  '.LLMDA,’  ’,FC, ’  ’, FCONSl,’  ' 

.PRI?rr*,’  ’,0,'  ', NOSEG 
RADCVT  -  PI/180.0 
DO  111  K10C-0,!!OS£G 
NUS ECO  -  NUSEG(KKK) 

NUSEG(KKK)  •  NUSEG(KKfC)  .*  RADCVT 
PRINT*,’  ’,NUSEG(KXK),’  ’ .MUSECD 
111  CONTINUE 
\^RITE(6,30) 

30  FORMATC ’ 1 ' , 7X , ' M ' , 8X , ' N ’ , 1 2X , 'REAL ’ , 14X , ’ IMACIOMARY ’ , 1 2Z , 

1  'MAGNITUDE ' , 14X , ' PHASE ' , 14X, 'TIME ' ) 

CALCL'LATE  CONSTANTS  FOR  USE  IN  THE  IIAIN  SUBROUTINE  AND  RELATED 
FUNCTIONS. 


FC0NS2  -  ((C**2)/2.0)  *  (COSII(2,0*HU)*SINIl(TAU)  +  TAU) 

AZERO  -  YNU(1.0S-03,0) 

AOI^E  -  -Yl(1.0S-03) 

DO  20  M-0, NOSEG 

DO  10  M  -  0, NOSEG 

CALL  CLZMNR(ZMNR,ZMNI) 

THE  RET'JR?rED  VALUES  FROM  CLZMf!R  ARE  REVERSED  DUE  TO  THE  PRE¬ 
SENCE  OF  THE  "J"  IN  THE  EXPRESSION  OF  THE  KERNAL.  THIS  IS 
L.VMDA  FOR  THE  FREDHOLM  INTEGRAL  EQUATION  OF  THE  SECOND  KIND. 


ZMN(M,N)  -  CMPLX(ZMNI,ZMNR) 

RMAC  -  SQRT(ZMMI**2  +  ZMNR**2) 

PHASE  -  ATAN(ZMNR/ZMNI)  *  180.0/PI 
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IF  (M.EQ.IJ)  THEN 

liTRITE (6,40)  M , N ,  ZMI.’I ,  ZMMR , RMAC ,  PHASE 
ELSE 

TOITE  (  6 , 50  )M ,  N,  Z;i?:i ,  ZML*R,  RrTAG,  PHASE 
END  IF 

40  F0?J1AT('0’,6(’*'),I3,6X,I3,4(6X,E15.8)) 
50  FORIIATC  ',6X,I3,6X,I3,4(6X,E15.8)) 

10  COOTINUE 

CALL  VOLTS (VMK) 

20  cormiRiE 

Ci\LL  CGECO(Z!frJ, MOSEG.NOSEG, IPVT.RCOjn), Z) 

PRItrr  *,’RCOND  -  '.RCOKD 

CALL  CGE  SL  ( ZMK ,  NOS  EC ,  NOS  EC ,  IPVT ,  V;i: : ,  0  ) 

WRITE  (8)  VIIN 

STOP 

END 


SL'BROUTINE  CLZMN?.(ZMKR,  ZMia) 

THIS  SUSnOLTINE  DOES  THE  CALCULATION  OF  TiiE  REAL  AND  TT!E  IMAG- 
IflARY  PARTS  OF  ZMH.  SINCE  THE  IMSL  INTEGRATOR  DOES  ONLY  REAL 
ARIT121ATIC,  THE  PARTS  HAD  TO  35  SEPERATED.  HENCE  THERE  IS  NO 
lIAIfKEL  FUNCTION  SL’BROLTINE  AS  HICIT  BE  EXPECTED. 

CO!  fMON/ELLIPS  /  IPJ ,  C ,  TAU ,  K 

COMMON/ CELLS /  !TU5 , NU6 , KUl , KU2 , NU3 , NU4 ,  ?  t , N, NOS  EG , NUSEG 
COMMON/COKST/IFLAC, FC, FCONSl , FI , PI .LAMDA 
EXTERiaiL  FAl,  FA3,  FA5,  FA7,  FAII,  FA3I,  FA5I,  FA7I 
EXTERNAL  FSIA 1 , FSIA2 ,  FSIA3, FSIA4 , FSI3 1 ,  FSIP. 2 , FS IB3 ,  FSIB4 
REAL  LAMDA , K, flU,  IRJSEC (0 ;  32)  ,NU1 , NU2 , NU3 ,  f;U4 ,  ;;U5 , NUS 


IFLAG  -  0  -  REAL  FUNCTION 

IFLAG  -  OTHER  -  IMAGINARY  FUNCTION 


PI2  -  PI/2.0 

FI  -  ((C**2)/2.0)  *  (C0SU(2.0*MU)*SINH(TAU)  +  TAU) 

IF  (M.GT.O)  THEN 

:iUl-MUSEG(M-l) 

ELSE 

NUl  -  -PI2  -  (PI2+NUSEG(l)) 

END  IF 
NU2-!IUSEC(M) 

IF  (M.EO.NOSF.C)  THEN 

NU3  -  PI2  +  (P12-NL'SEC(N0SEG-1)) 
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ELSE 

L'U3-MUSCC(M+1) 

END.  IF 

IF  (N.GT.O)  THEN 
i:U4-KUSEG(N-l ) 

ELSE 

:UI4  -  -PI2  -  (P12+NUSEG(1)) 

END  IF 
;:u5<:l’seg(n) 

I?  (N.LT. NOSEC)  THEN 
NU6-NUSEC(N+1) 

ELSE 

NU6  -  PI2  +  (P12-NUSEC(f:OSEG-l)) 

END  IF 

.  ITEST  -  M  -  N 

IF  ( ITEST. EQ.O)  THEN 
IFLAG  -  0 

ZM*rR  -(SriSON(FAl, NU2,NU3,NU5,NU6, 0.01, ERR, lER)  +  SD'S0N(FA5, 
IN  U1 , NU2 , NU4 ,  NT  5 , 0 . 0 1 , ERR , lER) 

\  +  3rtSON(FA3,*T2,tn]3,NU4,NU5,0.01,ERR,IER) 

A  +  SIMS0N(FA7,NU1,NU2,KU5,KU6, 0.01, ERR, ISR) 

1  +  SLMS0N(FA1I,NU2,NU3,PI-NU5,PI- 

2  0.01, ERR, lER)  +  SIMSOK(FA3I,NU2, 

3  NU3 , ri-NU4 , PI-NU5 , 0 . 01 ,ERR , lER 

4)  +  SI!1S0N(FA5I,NT1,NU2,P1-NU4,PI- 

5  MU5, 0.01, ERR, lEP)  +  SIMS0N(FA7I,NU1, 

5  NU2 , PI-NU5 , PI-NC6 ,0.01, ERR, 

7  lER))  *  LAMDA 

IFLAG-l 

ZMNI  -  Si:tSON{FAl,NU2,NU3,NU5,NU6, 0.01, ERR, lER) 

L  +  SIMSON(FA5,KU1,NU2,NU4,NU5,0.01,EP.R,IER) 

V  +  SI?IS0N(FA3, NU2,KU3,-NU4,NU5, 0.01, ERR,  lER) 

1  +  SIf’SON(FA7,KUl,HU2,NU5,KU6, 0.01, ERR,  lER) 

L  +  SirtSON(FAlI,rT2,NU3,PI-KU5, 

2  PI-NU6,0.1,ERR,ir.R)  +  SIMSON(FA3l, 

J  :RJ2  ,  NU3 ,  PI  -NU4 ,  PI-NU5 ,0.1, 

i  ERR.IER)  +  3I!lSOfI(FA5I,KUl,KU2,PI 

)  -NU4,P1-NU5,0.1,ERR,IER) 

>  +  SIHS0N(FA7I,NU1,NU2,PI-NU5, 

f  PI-KU6,0.1,ERR,IER) 

ZMKI  -  ZMNI  *  ( -LAMDA) 

L  +  FI  *  DCADRE(FSIA1,NU1,NU2,0.0,1.E-1,ERR,IER) 

!  -  TAU*C*C*DCA0RC(FSIR1,KU1,NU2,0.0,1.E-1,ERR,IER) 

>  .  +  FI  *  DCADRE(FSIA4,NU2,NU3,0.0,1.E-1,ERR,IER) 

>  -  TAU*C*C*DCADRE(FSIB4,KU2,NU3,0.0,I.E-1,ERR,IER) 


ELSE  IF  (ITEST.EO.l)  THEN 
IFLAG  -  0 

ZMNR  -  (SLMSON(FA7, NUl,NU2,NU5,r:U6, 0.01, ERR, IF.R)  +  SIMS0N(FA1, 
NU2,KU3, 
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1  KU5,i;U6, 0.01,  ERR,  lER) 

+  S  LM  SON  ( F  All ,  IR:2  ,  KL’3 ,  Pl-NU  5 , 

PI-riUf.,0.ni,ERR,IER)  +  5IMS0:I(FA3I, 

NU2 , NU3 , PI-NU4 ,?I-NU5, 

Q.01,ERr.,lER)  +  SriSON(FA3,NU2,NU3, 

J?U4,NU5, 0.01,  ERR.  lER)  +  SD:S0N(FA5, 
NU1,MU2,WJ4,KIJ5, 0.01,  ERR, 
lER)  +  SIMSOK(FA5I,NUl,f:U2,PI-NU4,PI-NU5, 

0.01, ERR, lER)  +  5IMS0N(FA71,N'J1,NU2, 

PI-NU5, PI-NU6, 0.01, ERR, lER))  *  LAMDA 
IFLAC=1 

ZMNI  -  SIMS0N(FA7,t:Ul,NU2,NU5,>lU6,0.01,IER,ERR)  + 

SIMS0N(  FAl ,  NU2 ,  f:U3 ,  IHJS , 

NU6,0.1,ERP.,IER)  +  SIMS0N(FA11, 
tn;2 ,  NU3 ,  PI-NU5 ,  PI-NU6 , 

0.1,ERR,ISR)  +  SIMSOi;(rA3,NL’2,NU3, 
NU4,NU5,0.1,EP.R,IER)  +  Sl.’IS0r;(FA3I , 

.•;U2 ,  !.nj3 ,  PI  -NU4 ,  PI-NU5 , 0 . 1 ,  ERR , 

IF.R)  +  SIHSON(FA5,KU1,KU2,j:U4,NU5,0.1,ERP., 
lER)  +  SIMSON(FA5I,t;Ul,NU2,PI-:iU4,PI-NU5,0.1, 

■  ERR.IER)  +  SIMS0N(FA7:I,>;U1,KU2,PI-Nn5, 
?I-KU6,0.1,ERR,IER) 

2MNI  -  ZMNI  *  (-LA.'tDA) 

+  Fl*DCADRE(FSIA2,:JUl,t.’U2,0.n,l.E-l, ERR.IER) 

-  TAU*C*C*DCADRE(  F3IB2 ,  NUl  ,;:U2 , 0 . 0 , 1 .  E-1 ,  ERR,  lER) 

ELSE  IF(ITEST.EQ.-l)  THEN 
IFLAC  -  0 

LMNR  -  (SIMS0K(FA3,NU2,NU3,NU4,in;5, 0.01, ERR.IER)  + 

S  lit  SOfJ(  FAl ,  NU2 ,  NU3 ,  IIU5 ,  NU6 , 0 . 01 ,  ERR 

,IER)  +  SIMSON(FAII,NU2,MU3,P1-NU5,PI-MU6,0.01, 
ERR.IER)  +  SIMS0N(FA3I,NU2,NU3,PI-NU4,PI- 
tRJ5, 0.01, ERR.IER)  +  Sl:iS0i:(FA5,NUl  ,NU2,:ru4,r;u5, 
0.01,  ERR.IER)  +  SLMSOK(FA5l,!:Ul,:nj2, 
PI-NU4,PI-NU5, 0.01, ERR.IER)  +  5IMSOt:(FA7, 
MUl,MU2,::U5,i;U6, 0.01,  ERR.IER)  +  SIMS0N(FA7I  ,NU1 , 
NU2,PI-NU5,PI-NU6, 0.01, ERR.IER))  *  LAMDA 
IFLAG  -  1 

:fTr;i  -  SIMSON(FA3,NU2,NU3,KU4,!IU5,0.01,ERR,IER) 

+  S  LMSO.'K  FAl ,  KU2 ,  i;U3 ,  KU5 ,  NU6 , 0 . 0 1 ,  ERR 

,IER)  +  SlriSO:;(FAlI,NU2,NU3,PI-NU5,PI-NU6,0.1, 
ERR.IER)  +  SIMSOri(FA3I,NU2,NL'3,PI-:iU4,PI- 
KU5, 0.1,  ERR.IER)  +  SIMSON(FA5,NUl,fIU2,NU4,rRJ5, 
0.1, ERR.IER)  +  SIMS0N(FA5I,NU1,NU2,PI-NU4, 

PI -NU5, 0.1, ERR.IER)  +  SIMS0N(FA7,NU1 ,NU2,NU5, 
NU6, 0.1, ERR.IER)  SIMS0N(FA7I,;RJ1  ,NU2,PI-NU5, 
PI-NU6, 0.1, ERR.IER) 

2Mi:i  -  ZMNI  *  (-LArtDA) 

1  +  Fl*nCADRE(FSIA3,NU2,NU3, 0.0,1. E-1, ERR.IER) 

2  -  TAU*C*C*DCADRE(FSIB3,NU2,NU3, 0.0,1. E-1, ERR.IER) 

no 


o  o  n  o 


.\pp&;n)ix  D 


MATRIX  CALCULATOR  AMD  LINEAR  ALGEBRA  EQUATION  SOLVER  PACE  D-  5 


CALCULATES  THE  DISTANCE  BETCEEN  TOO  POINTS  OF  COKSTAfT  MU 
OK  THE  ELLIPTIC  SHELL 
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REAL  PUMCTION  FI3:RXL(V,VP) 

THIS  FUKCTION  DETERMItXS  WHICH  (REAL  OR  IMACIKARY)  SECTIONS  IS 
BEING  CALCULATED  3Y  THE  DEFINITION  OF  IFLAC  SET  BY  THE  CALLING 
PROGRAM.  SEE  APPENDIX  C. 

com:  *0M/CELLS  /  NU5  ,  NU6  ,  NUl ,  NU2 ,  NU3 ,  im ,  M ,  N ,  NOSEC ,  NUS  EG 

COMMO!:/ELLIPS  /MU ,  C ,  TAU ,  K 

COMMON /CONST/ IFLAG , FC, FCONSl , FI ,P1 ,LAHDA 

COflHOK/S  INGLE  /  AXERO ,  AOtJE 

REAL  K ,  MU ,  rrUS EC  ( 0 : 32  ) , NUl , NU2 , NU3 ,  NU4 ,  NUS ,  NU6 

ARG  -  K  *  C  *  DEL(V,VP) 

IF  (IFLAG.CQ.O)  THEN 

FKERNL  -  (-1.0)  *  BESJ0(x\RC)  *  FKEPXKV.VP)  *  (C**2) 

ELSE 

IF  (ARG.LT.l.OE-03)  THEN 

SEED-  A2ER0  -  (AONE  *  5.0E-04> 

SEETNO  »  AO!'IE/2.0E-03 

FYJiRNL  -  -(SEED  +  SEE-n:0»(x\RG**2 ))  *  FXERL1(V,VP)  *  (C**2) 
ELSE 

FKERNL  •  -YNU(ARG,0)  *  FTSRLKV.V?)  *  (C**2) 

EN’d'  if  ' 

SIR)  I? 


RETURN 

ENT) 


REAL  FUNCTION  FI3:RL1(V,VP) 


COMMON/ELLIPS /MU, C , TAU, K 
COMMON/CONST/ IFLAG,  FC,  FCONS 1 ,  FI ,  PI ,  LA:>DA 
REAL  K,MU 

FUNCTION  CALCLOATES  THOSE  EXPRESSION’S  ASSOCIATED  l/ITH  THE  ZERO 
ORDER  IIANKEL  FUNCTIONS. 

FKERLl-(TAU**2)*(COSH(MU)**2-COS(V)**2)*(COSi:(MU)**2-COS(VP)**2) 


RETURN 


APPE;;DIX  D.  MATIIIK  calculator  and  Lirn^AR  ALCESr^A  equation  SOLV’ER  pace  d 


REAL  FUNCTION  BASEl(V,K,\Rll.‘:US,V>!) 

CALCULATES  THE  BASIS  FUNCTION  FOR  THE  VSUBMMINUSONE  TO  THE 
VSUB’I  TERM. 

COMMON/ELLIPS/MU.C 
REAL  K,  MU 

DIFWM  -  DEL(V,VMINUS) 

DIFVMV  -  DEL(\^t,VMINUS) 

BASEl  -  SIK(K*C*DIFVVM)  /  SIN(K*C*DIF\RIV) 

RETURN 

END 


RE,AL  FUNCTION  EASE2(V,K,VPLUS  ,V!!) 

FUNCTION  FOR  THE  PALLING  PORTION  OF  THE  SI>50SI[!AL  PORTION 
OF  TIP-  BASIS  CURVE  -  ItriEGRATE  FRai  VSUBM  TO  VSUSMPL'JSONE. 

COMMCNVELLIPS  /  MU ,  C ,  TAU 
REAL  K,MU 

D1FV\T  -  DCL<VPLUS,V) 

DIFVPV  -  DEL(VPLUS,VM) 

Bx\SS2  -  SIN(K  *  C  *  DIFVVP)  /  SIN(K  *  C  *  DIFVPV) 


RETURN 
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REAL  FUNCTION  FA3(VP) 

CO:iMOi:/CELLS/NU5 ,  N'Uf ,  NUl  ,NU2  ,NU3  ,i:U4 ,  !1,  N,  r:OSEG ,  rX'SEG 

CO!  IMON/ELLIPS  /  ftU ,  C ,  TAU ,  K 

CaiMON/VARY/V 

REAL  K, !fU, tX'SEC(0 : 32 )  ,NU1 , NU2, KU3 , r;U4 , NU5 ,  NU6 

FA3  -  BASEl(VP,K,NRl4,:nJ5)  *  BASE2(V,K,NU3, 

1  KU2)  *  FF3RrX(V,VP) 


RETURN 

END 


RE.\L  FUrcTTON  ROIISCE(V) 

CO’  fMON’/ELLIPS  /fRT ,  C ,  TAU ,  K 

cav.  !ot’/son?Ci:/  cs ,  iqiot,  vs 

RE.\L  K.NU.SCJOT 

ROly’SCS  -  SORT<(COSH(!rJ)**2)  *  ((C  *  COS(V)  ~  CS  *  COSD(VS)  )**2) 
I  +  (SINU(JJU)**2)  *  ((C  *  SIN(V)  -  CS  *  SINU(VS))**2)) 


RETL'RN 

END 


SUBROUTINE  VOLTS  (VMN) 

COHMON/CELLS/NU5 , NU6 , NUl , NU2 , KU3 ,IRJ4 ,M , N, NOSEC , NUSEC 
COMMON/CONST/IFLAC , FC , FGONS 1 , FI , PI , LAMDA 
REAL  KNOT,nUSEC(0:32) ,NUl,KU2,NU3,NU4,tIU5,NU6 
COMPLEX  VMN(0:32) 

EXTERNAL  PS4,FS5 

IFLAC  -  0 

VMNR  -  0CADRn(FS4,NUl,NU2,0.0,1.0E-3,ERR,IER 
1  )  +DCADRE(FS5,KU2,NU3,0.0,1.0E-3,ERR,IER) 

IFLAG  -  1 

VMNI  -  DC.\DRE(FS5,MU2,fnJ3,0.0,1.0E-3,ERR,IER 
1  )+DCADRE(FS4,NUl,NU2,0.0,1.0E-3,SRR,IER) 

VMii(M)  -  CMPLX(v?r:ni,viiNi) 

RIIAC  -  SQRT(VM!CR**2  +  VMNI**2) 

PHASE  -  ATAN(VIINI/VriNR)  *  (180,0/PI)  _ _ 

\^RITE(6,10)  H,VJtI!R,VRINI,RMAC,PHASE  [Rtproduetd  (rem  J|k 

10  F0RMAT('0',6X,I3,9X,4(6X,E15.8))  *v»iUbl»  copy.  1|P 
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RETURN 

END 


REAL  FUMCTIOM-  FA7(VP) 


COf  U  lOM/CELLS  /KU5 ,  r.’U6 ,  MUl ,  :,’U2 , 1;U3 ,  NU4 ,  ,  E,  KOSEC ,  XUS  EG 

C0M!10K/ELLIP5/>RJ,  C,  TAU.K 

CO!tl!Ot;/VARY/V 

REAL  K ,  raT ,  NTS  EC  (0 : 32  )  ,  tXl ,  ERJE ,  MU3 ,  MU4 ,  MU5 ,  MUfi 
FA7  -  BASE1(V,K,NU1,MU2) 

1  *  RASE2(VP,K,MU6,fJU5)  *  FIXR::L(V,VP) 


RETL'Ri: 

EKD 


REAL  FUECTIOK  FS2(V) 

C0f1M0!.7ELLIPS/MU,C,TAU,K 
COM!  10X7  SOUP.CE/  CS ,  KIIOT ,  VS 
CC3MM0M/C0NST/IFLAC,  FC,FC0!JS1,  FI ,  PI  .LAJJUA 
REAL  K,KNOT,!tU 

IF  (IFLAG.EO.O)  THEN 

FS2  •  FCOrai  *  (C**2)  *  TAU  *  BESJO(KMOT*RW:SCE(V) ) 

I  *  (C0SH(!IU)**2  -  C0S(V)**2) 

ELSE 

FS2  -  FCOKSl  *  YX’U((KNOT*ROt'SCE(V)),0)  *  (C**2)  *  TAU 
I  *  (C0SH(M!T)**2  -  C0S(V)**2) 

END  IP 


RETURN 
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REAL  FUNCTION  FSA(V) 

CO!  IMON/  CELLS  /:;US ,  KUf ,  HU  1 ,  MU2 .  rX’3 ,  NU4 ,  *t ,  K ,  I.’OSEG ,  trUSEC 
CCXIMOK/SOUP  CE/  CS ,  KIIOT 

REAL  KK0T,NUSEG(0: 32)  ,NUl,:iU2,NU3,NU4,!;U5,t;U6 
FS4  -  BASE2(V,rjmT,KU3,NU2)  *  FS2(V) 


RETURN 

ENT) 


REAL  FUNCTION  FS5(V) 

COtIHONVCELLS  /NU5 ,  KU6 ,  KUl ,  NU2 ,  t:U3 ,  !IU4 ,  !I ,  I! ,  KOSEC ,  NUS  EG 
COfl!  !Ot:/SOURCE/CS ,  lO.'OT 

REAL  KJ!0T,NL’SEG(0;32)  ,NU1,1RJ2,NX’3,NU4,HU5,NU6 

FS5  -  BASEl(V,K:iOT,f;Ul,!X'2)  *  FS2(V) 

RETL’RtJ 

EtU) 


REAL  FUNCTION  FAl(VP) 

CO!  II  ION/  CELLS  /  NUS ,  NU6 ,  NUl ,  NU2 ,  NU3 ,  NU4 ,  M ,  N ,  NOS  EG ,  NUS  EG 

COHMON/ELLIPS  /!IU ,  C ,  TAU ,  K 

CaiMON/VARY/V 

REAL  K,riU,NUSEC(0:32)  ,NU1,NX'2,NU3,NU4,NU5,NU6 

FAl  -  FKERNL(V,VP)  *  BASE2(V,K,NU3,NU2) 

1  *  BASE2(VP,r,,NU6,MU5) 


RETLT^N 
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REAL  FUNCTION  FAS (VP) 

C 

cw t:  ion/  cells  /  nus  ,  r:u6 ,  nui  ,  f:u2 ,  i:u3 ,  nu4  ,  m  ,  n  ,  nos  eg  ,  nus  eg 

COMMON/ELLirS / MU , C , TAU , K 
CaiMON/VARY/V 

REAL  K, IPJ,  tR;SEG(0 : 32 )  ,KU1 ,  NU2 , KU3 ,  KU4 , NUS ,  NU6 
C 

FAS  -  BASE1(V,K,MU1,NU2)  *  FKERBL(V,VP) 

1  *  BASEl(VP,K,NU4,f;US) 

C 

RETURN 

END 


REAL  FUNCTION  SIMSON(FUNSON,A,B,C,D,E,ErJlOR,IER) 

C 

EinERNAL  FUNS  ON 
<.:OfIMOK/VARY/Vl 
UF-AL  H,IJrr,FUNSOM 
C 

C  THIS  FL'NCTION  IS  BASED  ON  ROUTIL’E  GIVEN  IN; 

C 

C  'APPLIED  NUMERICAL  METHODS  FOR  DIGITAL  CaiPlTATION 

C  UITU  FORTRAN  AND  CSMP' , 'SECOND  EDITION',  BY  M.  L.  JAMES, 

C  G.M.  SMITH,  AMD  J.C.  UOLFORD,  NSW  YORK:  THOMAS  Y.  CROWELL. 

C  FACES  32C, 331, (1977). 

C 

n-(B-A)/lO. 

SUM-O.O 
■  Vl-A+a 
DO  10  I  -  2,10 
IF  (M0D(I,2))  20,20,30 
20  CONTINUE 

SUM-SU:i44  .*DCADRE(FLMISON,C,  D,  0.0, 1 .  E-1 ,  ERROR,  lER) 

GO  TO  10 
30  co?mtnjE 

SLPl  -  SUM-f2.*DCADRE(FUNS0!;,C,D, 0.0,1. E-1, ERROR, lER) 

10  VI  -  VI  +  H 

VI  ■  \ 

INT  -■  (11/3. 0)*(DCADRE(FUNSON,C,D,0. 0,1. E-1, ERROR, lER) 

1  +  SUM) 

VI  -  B 

INT  -  INT  +  (11/3.0)  *  DCAORE(FUKSON,C, D, 0.0,1. E-1, ERROR, lER) 
DELT  -  DEL(A,B) 

DELT2  -  0EL(C,D) 

C 

SIMSON  -  INT 
RETURN 


A?PS!n5i::  n. 
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END 


REAL  FUKCTIOII  FAII(VP) 

C0M?!0:,7CELLS  /MU5 ,  MU6 ,  IJUl » NU2 ,  KU3 ,  NU4 ,  f ! ,  K,  NOSEG ,  FUSEG 

CafMOW/ELLIPS/MU, C, TAU, K 

C0M^f0^7VAR\7V 

REAL  y. ,  MU,  IIUSEC (0 : 32  ) . MUl , KU2 , KU3 ,  KUA ,  NU5 ,  ::U6 

FAll  -  FKERIJL(V,VP)  *  BASE2(V,K,OT3,N*U2) 

1  *  BASEl(VP.K.^nJ6,I•IU5) 

RETCRN 

END 


REAL  FUNCTION  FA3I(VP) 

CaiMOrVCELLS  /N  U5 ,  ^6 ,  NUl ,  KU2 ,  NU3 ,  KUA , !  r ,  N ,  KOS  EC ,  KUS  EC 

CO!  !!10r:/  ELLIPS  /MU ,  C ,  TAU ,  K 

COJnOfl/VARY/V 

REAL  K , !  lU ,  NUS  EC  <  0 : 3  2  ) ,  NUl ,  KU2 ,  NTJ3 .  MUA ,  NU  5 ,  ::U6 

FA31  -  FKEP.NL(V,VP)  *  BASE2(V,K,NU3,NU2) 

1  *  BASE2(VP,K,NUA,NU5) 


C 


RETURN 

END 


k 
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REAL  FUNCTION  FA5I(VP) 

coriMON/CELLS  /NU5 ,  r:u6 ,  i:ui .  :pj2  ,  i:u3 ,  ku4  ,  ri , ,  nosec  ,  uus  eg 

COM!  lOI!  /ELLIPS  /.MU ,  C ,  TAU ,  K 

CO:itlOK/VARY/V 

REAL  K ,  !1U ,  NUS EG ( 0 : 32  ) , NTIl , NU2 ,  rnJ3 , NUA ,  NU5 ,  KU6 

FA5I  -  FKERKL(V,VP)  *  BASE1(V,K,KU1 ,KU2) 

1  *  EASE2(VP,K,NU4,f.'U5) 

RETURN 

END 


REAL -FUNCTION  F.A7I{VP) 

CCrjMOM/CELLS  /IJU5 ,  r7U6 ,  IMIl ,  rnJ2 ,  NU3 ,  -NUA ,  M ,  N,  KOSEG ,  NUS  EG 

COMMOt:/ELLIPS  /MU ,  C ,  TAU,K 

C0!IM0I,7VARY/V 

REAL  K,.*IU,NUSEG(0; 32) ,KUl,NU2,!nJ3,riUA,KU5,nU6 

FA7I  -  FKERIIL(V,VP)  *  RASEI(V,X.NU1,KU2‘) 

1  *  BASE1(VP,K,KU6,NU5) 

RETURN 

END 


REAL  FUNCTION  FSIAl(V) 

CC3H!  10N/CELLS  /t!U5 ,  ::U6 ,  NUl ,  NU2 ,  NU3 ,  CTA ,  M ,  t! ,  !:OS  EG ,  NUS  EG 
COMMON/ELLIPS  /?1U ,  C ,  TAU  ,  K 

REAL  K ,  MU ,  NUS  EC  (0 ;  32  )  .l.'Ul ,  NU2 ,  NU3 ,  NUA ,  MU5 ,  f  :U6 

FSIAl  -  BASE1(V,R,NI!1,:TU2) 

1  *  EASE1(V,K,NU4,NU5) 

RETURN 

END 
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RE.\L  FUNCTION  FSIBl(V) 


COMMON/CELLS  /  MU5 ,  ?:Ufi ,  KUl ,  NU2 ,  KU3 ,  NU4 ,  N ,  N ,  NOS  EC ,  NUS  EG 
CO?  1M0N7ELLIPS  /  .'fU ,  C ,  TAU ,  K 

REAL  K,MU,NUSEG(0;32),!nJl,NU2,i:U3,NU4,NU5,N’U6 

FSIRl  -  FSIAl(V)  *  (COS(V)**2) 

RETURN 

E?ro 


REAL  FUNCTION  FSIA2(V) 

COM?^Ot:/CELLS  /i:U5 ,  JX'6 ,  NUl ,  NU2 ,  fnJ3 ,  NU4 ,  M ,  N,  frOSEC ,  NUS  EG 
CO:  IMON/ELLIFS  /'?U ,  C ,  TAU ,  K 

REAL  K ,  MU ,  NUS  EC  (  0 ;  32  ) ,  NUl ,  KU2 ,  NU3 ,  NU4 ,  NUS ,  i:U6 

FSIA2  -  EASEl(V,i;,m,MU2) 

1  *  BASE2(V,K,NU6,:JU5) 

RETURN 

END 


REAL  FUNCTION  FSIB2(V) 

COtI‘fOf;/CELLS  /NU5 ,  NU6 ,  NUl ,  NU2 ,  KU3 ,  NU4 ,  M ,  N,  NOSEC ,  NUS  EC 
CO;iHON/ELLIPS/rfU,C,TAU,K 

REAL  tIU,K,NUSEC(0;32)  ,NUl,rRJ2,NU3,NU4,!a:5,NU6 

FSIB2  -  FSIA2(V)  *  (COS(V)**2) 

RETURN 

END 
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REAL  FUNCTION  FSL\3(V) 

COMMOKVCELI^  /NU5 ,  NU6 ,  UUl ,  KU2 ,  NU3 ,  KUA ,  U,  KOSEC,  NUSEG 

con?!or;/ELLirs/trj,  c,  tau,  k 

REAL  MU , K, NUS EG ( 0 : 32 ) , NUl ,  ND2 ,  «D3 , MUA , KU5 , NU6 

FSIA3  -  3ASE2(V,K,Nl'3,KU2) 

I  *  BAS£1(V,K,KUA,NU5) 

REILDN 

END 


REAL  FUlXTION  FSIB3(V) 

COMMO!!/CELLS  /NU5 ,  NU6 ,  NUl ,  ?n:2 ,1X3 ,  IXA ,  M,  N ,  KOSEG ,  NUS  EG 
COM!  lOM/ELLIPS  /IX ,  C ,  TAU ,  K 

REAL  MU , K, NUSEG (0:32)  ,NU1 , 1X2 , :1U3, NUA , MU5 ,  NX6 

FSIB3  -  FSIA3(V)  *  CnS(V)**2 

RETURN 

END 


REAL  FUNCTION  FSIAA(V) 

C 

COtlMOK/CELLS  /KU5 , 1X6 ,  NUl ,  1X2 , 1X3 ,  NUA  ,  M,  N,  NOS  EG ,  NUS  EG 

com:  ion/ellips  /ix  ,  c  ,  tau  ,  k 

REAL  MU, K,  1XSEC(0 : 32 ) ,NUl ,NU2 , NU3, NUA , NUS , NU6 

FSUA  -  BASE2(V,K,NU3,NU2) 

1  *  BASE2(V,K,N'U6,KU5) 

C 

RETURN 

Eiro 


C 

C 

C 
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REAL  FUNCTION  FSIR4(V) 

CO!!MOi:/CELLS /NU5 ,  NU6 ,  NUl , NU2  ,  t:U3 , NU4 ,  M ,  N,  NOSEC ,  NUSEC 
COIIMON/ELLIPS  /MU ,  C,  TAU ,  N 

REAL  MU ,  K ,  NUS EG (0:32),  NUl ,  MU2 ,  KT3 ,  NU4 , JOJS ,  MU6 
FSIB4  -  FSIA4(V)  *  C0S(V)**2 


RErJRN 


r--’ 


a 


c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 


c 

c 

c 

c 


REAL  FU::CTtO:i  Yl(ARG) 

SITRROLTINC  FRO»>  'SYSTEM/350  SCIENTIFIC  SUBROlTItni  PACKAGE 

(36nA-cri-03y.)  versio?:  ii,  procram’JErs  ma-mual',  .120-0205-2, 

I.'rrCRNATIONAL  BUSIN’ESS  MACHINES  (IBM)  CORPORATION’,  WHITE 
PLAINS,  ND./  YORK.  PAGES  157  &  15S,  1957. 

SEE  ALSO  PAGE  275  FOR  ACCUPvACY  IKFOP.MATION. 

C0M?!0N7C0HST/ I F LAG ,  FC ,  FCOnS  1 ,  FI ,  pi 
IF  (ARG.GE.4.)  THEN 
XI  -  4.0  /  ARG 

T2  -  16.0  /  ARG**2 

PI  .  ((((4.2414E-06*T2-2.00920E-05)*T2+5.80759E-05)*T2 

1  -  2.232030E-04)*T2+2.9218256E-03)*T2+0.3?S9422P19 

Ql  -  ((((-3.6594E-06*T2+1.62200E-05)*T2-3.9R708E-05)*T2 

1  +  1.064741E-04)*T2-6.390400E-04)*T2+3.7400S36E-02 

NOTE;  ERRORS  WERE  NOTED  IN  THE  IBM  WRITE  UP  ON  PACE  58.  CODE 
USED  IS  CORRECT.  SEE  HITCECOCK,  A.  J.  M.  -POLYtrOMIAL  APPROX¬ 
IMATION’S  TO  3ESSEL  FUNCTIONS  OF  ORDER  ZERO  AND  0^:E  AND  TO 
REL/XTED  FUNCTIONS".  7UTEMATICAL  TABLES  AND  OTHER  AIDS  TO 
CaiPUTATION ' ,  'VOLUME  ll’(58),  PAGES  86-88,  APRIL,  1957. 

.AA-2.0/SORT(ARG) 

BB-AA*T1 

D  -  ARG  -  PI/4.0 

Yl— AA*P1*C0S(D)+RB*Q1*SIN(D) 

ELSE 

XX-ARG/2.0 

X2*XX^XX 

*T  -ALOG(XX)+0. 5772156649015 

EULER'S  CONSTAtrr  FROM  'HANDBOOK  OF  TABLES  FOR  MATHEMATICS', 
'THIRD  EDITION',  PAGE  5,  1967. 


SUM  -  0.0 

i'J 

APPENDIX  D.  MATRIX  CALCLUVTOR 

« 

!•/ 

. 

TER.M  -  XX  *  (T  -  0.5) 

YONE  -  TER^f 

DO  80  L  -  2,16 

SUM  -  SU!I  +  1.0/  FLOAT(L- 

FL  -  FLOAT(L) 

%  • 

FLl  -  FL  -  1.00 

■|  * 

•  TS  -  T  -  SUM 

•*- 

TERJl  -  (TERM*(-X2)/(FL1*FL))* 

YON’S  -  YONE  +  TERM 

s  ^ 

80  CONTINUE 

PI2  -  2.0  /  PI 

Y1  ■  (-P12/AP.G)  +  PI2  *  YOr® 

END  IP 

RETURN 

! 

END 

»  ' 

« 

It'*. 

-  17 
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